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ABSTRACT 

The  High  Enthalpy  Shock  Tunnel  Gottingen  (HEG)  of  the  German  Aerospace  Center  (DLR)  is  one  of  the 
major  European  hypersonic  test  facilities.  It  was  commissioned  for  use  in  1991  and  was  utilised  since  then 
extensively  in  a  large  number  of  national  and  international  space  and  hypersonic  flight  projects. 
Originally,  the  facility  was  designed  for  the  investigation  of  the  influence  of  high  temperature  effects  such 
as  chemical  and  thermal  relaxation  on  the  aerothermodynamics  of  entry  or  re-entry  space  vehicles.  Over 
the  last  years  its  range  of  operating  conditions  was  subsequently  extended.  In  this  framework  the  main 
emphasis  was  to  generate  test  section  conditions  which  allow  investigating  the  flow  past  hypersonic  flight 
configuration  from  low  altitude  Mach  6  up  to  Mach  10  in  approximately  33  km  altitude.  The  studies 
performed  focused  on  the  external  as  well  as  internal  aerodynamics  including  combustion  of  hydrogen  in 
complete  supersonic  combustion  ramjet  configurations.  The  complexity  of  these  flows  requires  that 
experiments  in  ground  based  facilities  are  strongly  linked  with  computational  fluid  dynamics  (CFD) 
investigations.  These  common  activities  range  from  the  calibration  process  of  the  facility  and  the  study  of 
basic  aerodynamic  configurations,  which  are  well  suited  to  look  at  fundamental  aspects  of  high  enthalpy 
flow  fields,  to  the  investigation  of  complex  re-entry  configurations.  In  the  DLR  Spacecraft  Department  the 
research  programs  conducted  in  HEG  are  closely  linked  with  numerical  studies  using  the  DLR  TAU  code. 


1.0  INTRODUCTION 

In  hypersonics,  the  high  cost  and  risk  of  designing  complex  vehicles  precludes  extensive  experimental 
prototype  flight  testing.  Therefore,  ground  based  testing  facilities  were  developed  and  they  played  an 
important  role  since  the  early  era  of  hypersonic  flight.  Later,  as  computing  resources  became  more 
advanced,  computational  fluid  dynamics  (CFD)  tools  were  developed  and  utilised.  Recent  approaches 
undertaken  in  the  framework  of  projects  such  as  HyShot  [66]  or  SHEFEX  [52]  show  that  coupling  of  the 
three  main  tools  of  hypersonic  flight  vehicle  design,  namely  hypersonic  ground  based  testing,  CFD  and 
flight  testing,  are  becoming  economically  achievable  (Figure  1).  The  means  to  arrive  at  this  goal  is  to 
perform  hypersonic  flight  tests  of  new  technologies  based  on  sounding  rocket  technology  which  can  be 
realized  with  an  order  of  magnitude  lower  budget  compared  to  flying  complex  X- vehicles.  This 
development  will  lead  in  the  near  future  to  more  frequent  flight  tests. 

The  objectives  of  the  present  article  is  the  discussion  of  the  closely  coupled  experimental  and  numerical 
approach  for  hypersonic  and  high  enthalpy  flow  investigations  utilising  the  DLR  High  Enthalpy  Shock 
Tunnel  Gottingen,  HEG  and  the  DLR  TAU  Code.  While  flight  testing  is  the  ultimate  goal  and  required  to 
close  the  loop  shown  in  Figure  1,  it  is  only  briefly  addressed  here. 
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Numerical  analysis  (DLR  Tau  code) 


Figure  1 :  Three  major  sources  of  aerothermodynamic  data  -  ground  based  testing,  numerical 

analysis  and  flight  experiment 

The  article  is  structured  such  that  after  a  description  of  the  HEG  facility  (chapter  2)  and  the  applied 
measurement  techniques  (chapter  3),  the  DLR  TAU  code  is  presented  in  chapter  4.  Finally  in  chapter  5, 
examples  of  combined  experimental  and  numerical  flow  field  investigations  are  presented.  These  comprise 
the  calibration  of  the  HEG  free  stream  flow,  high  enthalpy  cylinder  shock  layer  investigations,  free  jet 
testing  and  CFD  analysis  of  the  HyShot  II  scramjet  flight  experiment  and  the  HEG  post  flight  analysis  of 
the  DLR  SHEFEX-I  flight  experiment. 


2.0  EXPERIMENTAL  TOOLS 

2.1  Hypersonic  and  High  Enthalpy  Ground  Based  Testing 

In  hypervelocity  flows  the  speed  of  the  considered  fluid  is  much  larger  than  the  speed  of  sound. 
Commonly  the  hypersonic  flow  regime  is  considered  to  start  above  a  Mach  number  of  M=5.  Ground  based 
testing  of  such  flows  is  performed  in  many  different  types  of  facilities.  The  reason  for  this  is  the  large 
range  of  flow  conditions  and  phenomena  encountered  in  hypersonic  flight  and  the  fact  that  no  single 
facility  can  simulate  all  relevant  flow  parameters  simultaneously.  Therefore,  in  hypersonic  testing,  partial 
simulation  of  different  flow  phenomena  is  performed  in  different  types  of  facilities.  Examples  are  Mach- 
Reynolds  number  simulations  in  cold  hypersonic  ground  based  test  facilities,  verification  and  qualification 
of  hot  structures  of  space  vehicles  in  arc -heated  test  facilities  or  the  investigation  of  the  influence  of  the 
chemically  reacting  flow  past  an  entry  or  re-entry  vehicle  on  its  aerodynamic  behaviour  in  shock  tunnels 
or  shock  expansion  tunnels.  Comprehensive  overviews  of  ground  based  testing  of  hypersonic  flows  are 
given  by  e.g.  Lukasiewicz  [50]  and  Lu  &  Marren  [49]. 

One  possibility  to  increase  the  Mach  number  in  ground  based  facilities  is  by  reducing  the  free  stream 
temperature,  i.e.  the  free  stream  speed  of  sound.  In  this  case,  high  Mach  numbers  can  be  achieved  while 
the  free  stream  velocity  is  significantly  lower  than  the  actual  flight  velocity.  However,  characteristic  of 
high  Mach  number  hypersonic  flight  with  M  -10  and  higher  is  that  the  kinetic  energy  of  the  flow  is  large 
enough  that  high  temperature  effects  such  as  vibrational  excitation  or  dissociation  of  the  fluid  molecules 
occur  in  the  flow  past  hypersonic  vehicles.  The  high  flow  velocities  and  subsequently  the  high  temperature 
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effects  are  not  duplicated  in  cold  hypersonic  ground  based  test  facilities.  During  the  re-entry  flight  of  a 
space  vehicle  in  the  earth’s  atmosphere  or  during  the  interplanetary  atmospheric  entry  of  space  vehicles  or 
meteorites,  speeds  in  excess  of  6  km/s  are  achieved.  Considering  a  flow  with  this  speed  in  a  wind  tunnel 
test  section  with  an  area  of  1  m2  and  a  density  of  0.003  kg/m3,  a  power  requirement  of  already  300  MW 
results.  Therefore,  continuous  flow  facilities  are  not  a  practical  way  to  generate  such  high  enthalpy, 
hypersonic  flows.  Additionally,  the  correct  simulation  in  ground  based  testing  of  the  chemical  relaxation 
length  of  the  dissociation  reactions  of  the  fluid  molecules  occurring  for  example  behind  the  strong  bow 
shock  in  front  of  the  nose  of  a  re-entry  vehicle,  requires  the  duplication  of  the  flight  binary  scaling 
parameter,  the  product  of  the  free  stream  density  p  times  a  characteristic  flow  length  L  (see  e.g.  Stalker 
[85]).  Consequently,  the  smaller  the  scale  of  the  wind  tunnel  model  is  chosen,  the  higher  the  free  stream 
density  or  pressure  needs  to  be.  Considering  the  flight  trajectory  range  of  a  re-entry  vehicle  in  about  70  km 
altitude,  where  typically  the  highest  heat  flux  occurs,  the  atmospheric  density  is  approximately  10"4  kg/m3. 
Using  a  geometrical  scaling  factor  of  30,  a  free  stream  density  in  the  ground  based  facility  of  0.003  kg/m3 
is  required.  If  a  flow  with  this  free  stream  density  and  a  velocity  of  6  km/s  is  generated  by  expansion  in  a 
convergent  /  divergent  hypersonic  nozzle  from  a  reservoir  at  rest  without  adding  energy,  a  total  specific 
enthalpy  of  about  23  MJ/kg  and  a  nozzle  reservoir  pressure  in  the  order  of  90  MPa  is  required.  This  results 
in  a  nozzle  reservoir  temperature  of  about  10000  K.  It  is  clear  that  such  conditions  can  only  be  achieved  in 
impulse  facilities  with  short  flow  duration.  The  most  successful  facility  types  which  are  able  to  generate 
high  enthalpy  and  high  pressure  hypersonic  flows  are  shock  tunnels  and  shock  expansion  tunnels  with 
typical  test  times  of  approximately  5  milliseconds  and  less.  It  should  be  emphasised  here  that  the  test  time 
scales  with  the  size  of  the  facility.  The  number  quoted  here  is  typical  for  facilities  with  shock  tube 
diameters  in  the  order  of  100  -  200  mm  and  nozzle  exit  diameters  of  less  than  approximately  1  m.  The 
principle  of  these  facilities  is  to  store  the  energy  over  a  long  period  of  time,  therefore  reducing  the 
necessary  power  requirement  and  subsequently  releasing  the  stored  energy  rapidly.  Due  to  the  high  flow 
speeds,  test  times  in  the  order  of  a  few  milliseconds  are  still  sufficient  for  the  development  of  a  steady 
flow  over  a  model.  According  to  Hornung  [38],  a  reasonable,  conservative  correlation  of  the  necessary  test 
time  to  establish  a  steady  flow  is  x  =  20  L/Uoo,  where  L  is  the  model  length  and  Uoo  is  the  free  stream 
velocity.  For  a  test  using  the  above  mentioned  flow  condition  and  a  0.3  m  long  wind  tunnel  model,  the 
required  test  time  would  be  1  ms.  The  high  pressure,  high  velocity  flows  which  can  be  generated  in  shock 
tunnels  and  shock  expansion  tunnels  makes  these  facilities  not  only  suitable  for  the  investigation  of  space 
vehicle  aerothermodynamics  but  also  for  studying  complete  airbreathing  propulsion  systems,  particularly 
supersonic  combustion  ramjets  (scramjets)  at  flight  Mach  numbers  of  M  =  8  and  above.  In  this  framework 
it  is  important  that  in  addition  to  the  free  stream  Mach  and  Reynolds  number,  the  correct  static  pressure 
and  temperature  are  established  at  the  combustor  inflow.  Further,  if  hypersonic  flight  configurations  are 
considered  which  can  be  tested  at  1 : 1  scale,  the  flight  free  stream  conditions  can  be  duplicated  in  these 
facilities  generating  the  same  pressure  and  similar  heat  flux  loads  as  experienced  in  flight.  In  the 
subsequent  sections,  the  operating  conditions  realised  in  the  High  Enthalpy  Shock  Tunnel  Gottingen,  HEG 
will  be  presented. 

2.2  High  Enthalpy  Shock  Tunnel  Gottingen,  HEG 

The  HEG  is  a  free  piston  driven  shock  tunnel  ([28],  [32])  which  was  developed  and  constructed  in  the 
framework  of  the  European  HERMES  program  over  the  period  1989  -  1991.  It  was  commissioned  for  use 
in  1991,  at  that  time  being  the  largest  facility  of  its  type  worldwide.  Since  then  it  was  extensively  used  in  a 
large  number  of  national  and  international  space  and  hypersonic  flight  projects.  The  research  activities 
which  were  always  strongly  linked  with  CFD  investigations  range  from  the  calibration  process  of  the 
facility  and  the  study  of  basic  aerodynamic  configurations,  which  are  well  suited  to  investigate 
fundamental  aspects  of  high  enthalpy  flows  to  the  investigation  of  complex  entry,  re-entry,  hypersonic 
flight  and  integrated  scramjet  configurations. 


RTO-EN-AVT-1 86 


8-3 


A  Closely  Coupled  Experimental  &  Numerical  Approach  for  Hypersonic  &  High 
Enthalpy  Flow  Investigations  Utilising  the  HEG  Shock  Tunnel  &  the  DLR  TAU  Code 


e-*a‘r±A-*irtv‘ 


In  a  free  piston  driven  shock  tunnel,  the  conventional  driver  of  a  shock  tunnel  is  replaced  by  a  free  piston 
driver.  This  concept  was  proposed  by  Stalker  [84].  A  schematic  and  wave  diagram  of  this  type  of  facility 
is  shown  in  Figure  2. 


Figure  2:  Schematic  and  wave  (x-t)  diagram  of  a  free  piston  driven  shock  tunnel 


Figure  3:  Schematic  of  the  High  Enthalpy  Shock  Tunnel  Gottingen,  HEG 


Free  piston  driven  shock  tunnels  consist  of  a  secondary  reservoir,  a  compression  tube,  separated  from  an 
adjoining  shock  tube  via  the  primary  diaphragm,  and  a  subsequent  nozzle,  test  section  and  dump  tank.  A 
schematic  and  photos  of  HEG  are  given  in  Figure  3  and  Figure  4,  respectively. 

The  high  pressure  air  stored  in  the  secondary  reservoir  is  utilised  to  accelerate  a  heavy  piston  down  the 
compression  tube.  During  this  quasi-adiabatic  compression  and  heating  of  the  light  driver  gas  (typically 
helium  or  a  helium  argon  mixture),  the  piston  reaches  a  maximum  velocity  in  the  order  of  300  m/s.  The 
driver  gas  temperature  increases  with  the  driver  gas  volumetric  compression  ratio.  When  the  main 
diaphragm  burst  pressure  is  reached  it  ruptures  and  the  wave  process  as  in  a  conventional  reflected  shock 
tunnel  is  initiated  (see  Figure  2).  A  shock  wave  is  moving  into  the  driven  section  and  the  head  of  a  centred 
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expansion  wave  is  moving  into  the  high  pressure  region.  The  numbers  used  in  Figure  2  denote  distinct 
regions  of  the  flow.  Region  1  contains  the  test  gas  at  the  initial  shock  tube  filling  conditions  and  region  4 
contains  the  hot,  compressed  driver  gas  after  piston  compression.  Region  2  contains  the  shock  compressed 
test  gas,  while  in  region  3,  the  driver  gas  processed  by  the  unsteady  expansion  wave  is  contained.  The  test 
and  driver  gas  are  separated  by  a  contact  surface. 


Figure  4:  Photographic  views  of  the  High  Enthalpy  Shock  Tunnel  Gottingen,  HEG 

After  reflection  of  the  incident  shock  wave  at  the  right  end  wall  of  the  shock  tube,  the  test  gas  is  brought  to 
rest  in  region  0.  Subsequently,  the  reflected  shock  wave  penetrates  the  contact  surface.  Depending  on  the 
local  conditions,  three  types  of  shock  wave  /  contact  surface  interaction  can  be  differentiated.  Due  to  the 
fact  that  the  shock  compressed  and  heated  slug  of  gas  in  region  0  is  used  in  reflected  shock  tunnel 
operation  as  the  reservoir  driving  the  flow  in  the  nozzle  and  test  section,  shock  tube  operation  in  tailored 
interface  mode  is  most  desirable  (Figure  5).  At  this  condition  the  pressure  in  region  0  remains  constant. 


press 

sure 

[\ 

press 

sure 

press 

sure 

Figure  5:  Shock  wave  /  contact  surface  interactions  for  undertailored  (left),  tailored  (middle)  and 
overtailored  (right)  interface  condition;  wave  (x-t)  diagram  (upper  row)  and  time  history  of 

pressure  in  region  0  (lower  row) 


For  undertailored  or  overtailored  interface  conditions,  the  pressure  in  region  0  is  decreasing  or  increasing, 
respectively  after  interaction  of  the  reflected  shock  with  the  contact  surface  (Figure  5).  Reflected  shock 
tunnels  are  characterised  by  a  convergent  -  divergent  nozzle  which  is  attached  to  the  end  of  the  shock  tube. 
A  thin  secondary  diaphragm  is  placed  at  the  nozzle  entrance  in  order  to  allow  evacuation  of  the  nozzle, 
test  section  and  dump  tank  before  the  run.  The  nozzle  entrance  diameter  is  chosen  sufficiently  small  such 
that  the  incident  shock  wave  is  almost  completely  reflected.  The  stagnant  slug  of  test  gas,  generated  by  the 
shock  reflection  in  region  0  is  subsequently  expanded  through  the  hypersonic  nozzle.  The  nozzle  flow 
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starting  process  is  characterised  by  a  wave  system  which  passes  through  the  nozzle  before  a  steady  flow  is 
established  (see  Figure  2).  The  incident  shock  wave  (a)  is  followed  by  a  contact  surface  (b),  an  upstream 
facing  secondary  shock  wave  (c)  and  the  upstream  head  of  an  unsteady  expansion  (d).  The  trajectory  of 
the  piston  is  chosen  in  a  way  that  after  main  diaphragm  rupture,  the  pressure  and  temperature  of  the  driver 
gas  in  region  4  is  maintained  approximately  constant.  This  is  achieved  by  selecting  the  velocity  of  the 
piston  at  diaphragm  rupture,  and  therefore  the  subsequent  movement  of  the  piston  such  that  it  compensates 
for  the  loss  of  the  driver  gas  flowing  into  the  shock  tube.  For  that  reason,  in  contrast  to  the  constant 
volume  driver  of  conventional  shock  tunnels,  the  free  piston  driver  is  a  constant  pressure  driver.  Due  to  the 
large  forces  occurring  during  the  operation  of  the  free  piston  driver,  the  compression  tube,  shock  tube, 
nozzle  assembly  is  allowed  to  move  freely  in  axial  direction.  An  inert  mass  placed  at  the  compression  tube 
/  shock  tube  junction  can  significantly  reduce  the  recoil  motion  of  the  facility  during  operation.  The  test 
section  and  the  dump  tank  remain  stationary.  A  sliding  seal  is  used  at  the  nozzle  /  test  section  interface. 

The  overall  length  of  HEG  is  62  m  and  it  weighs  280  t.  Approximately  a  third  of  its  weight  is  contributed 
by  the  inert  mass  (see  Figure  3  and  left  picture  of  Figure  4).  The  compression  tube  is  closed  by  a  hydraulic 
oil  system  (quick  disk  connect)  at  the  main  diaphragm  station.  The  shock  tube  is  connected  to  the  nozzle 
of  the  tunnel  at  the  downstream  closure,  which  is  also  driven  by  oil  hydraulics  to  close  and  seal  the  tunnel. 
The  compression  tube  has  a  length  of  33  m  and  a  diameter  of  0.55  m.  The  shock  tube  is  17  m  long  with  a 
diameter  of  0.15  m.  The  HEG  was  designed  to  provide  a  pulse  of  gas  to  a  hypersonic  convergent  - 
divergent  nozzle  at  stagnation  pressures  of  up  to  200  MPa,  and  stagnation  enthalpies  of  up  to  23  MJ/kg. 
Regarding  the  test  gas,  no  basic  limitations  exist.  The  operating  conditions  presented  in  the  present  article 
are  related  to  the  test  gas  air.  Additionally,  operating  conditions  using  nitrogen  and  carbon  dioxide  exist. 

2.2.1  HEG  Operating  Conditions 

Originally,  HEG  was  designed  for  the  investigation  of  the  influence  of  high  temperature  effects  such  as 
chemical  and  thermal  relaxation  on  the  aerothermodynamics  of  entry  or  re-entry  space  vehicles.  As 
discussed  above,  in  order  to  correctly  model  the  chemical  relaxation  occurring  behind  the  bow  shock  of  a 
re-entry  vehicle,  the  flight  binary  scaling  parameter  must  be  reproduced  in  ground  based  testing.  Further, 
for  high  enthalpy  testing  an  additional  driving  parameter  which  must  be  reproduced  is  the  flow  velocity. 
Therefore,  the  operating  conditions  of  HEG  are  first  discussed  in  Figure  6  in  terms  of  the  binary  scaling 
parameter  pL  and  the  flow  velocity  u.  Here  L  represents  the  length  of  the  considered  configurations.  In 
addition  to  the  HEG  operating  conditions,  the  most  important  fluid  mechanical  and  chemical  processes 
occurring  during  re-entry  of  a  spacecraft  in  the  Earth’s  atmosphere  are  depicted  in  Figure  6.  Further,  as 
reference,  the  flight  trajectories  of  a  lifting  body  re-entry  from  low  Earth  orbit  (IX V),  a  ballistic  super 
orbital  re-entry  (Apollo  11)  and  a  hypersonic  flight  experiment  (SHEFEX)  are  provided.  An  indication  of 
the  corresponding  flight  altitudes  is  given  in  the  right  diagram  of  Figure  6  showing  the  temperature 
variation  of  the  Earth’s  atmosphere.  The  transitions  between  regimes  of  different  physical  and  chemical 
properties  shown  in  Figure  6  depend  on  the  chosen  reference  length  and  vary  when  different 
configurations  are  considered.  Further,  the  boundaries  shown  have  only  symbolic  character.  In  reality,  no 
clear-cut  dividing  lines  between  the  different  regimes  exist.  The  Knudsen  number  given  in  Figure  6  shows 
that  the  HEG  operating  conditions  are  located  in  the  continuum  flow  regime.  The  high  energy  content  of 
re-entry  flows  leads  to  strong  heating  of  the  air  in  the  vicinity  of  a  spacecraft.  Depending  on  the 
temperature  level  behind  the  shock  wave  (i.e.  the  flight  velocity),  the  vibrational  degrees  of  freedom  of  the 
air  molecules  are  excited  and  dissociation  reactions  of  oxygen-  and  nitrogen  molecules  may  occur. 
Further,  ionisation  of  the  air  constituents  occurs.  The  high  temperature  effects  described  here  are  enabled 
by  energy  transfer  from  the  translational  energy  stored  in  the  random  motion  of  the  air  particles,  which  is 
increased  by  the  gas  heating,  to  other  forms  of  energy.  Because  this  energy  transfer  is  realised  by  air 
particle  collisions,  it  requires  a  certain  time  period  to  develop.  The  time  required  to  reach  an  equilibrium 
condition,  is  e.g.  defined  by  the  local  temperature  and  density.  Therefore,  depending  on  the  ratio  of  the 
relaxation  time  to  a  characteristic  timescale  of  the  flow,  the  chemical  and  thermal  relaxation  processes  can 
be  either  in  non-equilibrium  or  in  equilibrium.  Further,  along  a  re-entry  trajectory,  the  Reynolds  number 
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varies  over  several  orders  of  magnitude.  In  high  altitude  flight  the  wall  boundary  layer  of  a  re-entry 
vehicle  is  initially  laminar.  After  exceeding  a  critical  Reynolds  number  (shown  exemplarily  for  the  IXV 
configuration  in  Figure  6),  the  transition  from  a  laminar  to  a  turbulent  boundary  layer  takes  place.  This 
process  is  linked  with  an  increase  of  the  skin  friction  and  the  wall  heat  flux. 
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Figure  6:  HEG  operating  conditions  in  terms  of  the  binary  scaling  parameter  pL  and  the  flow 

velocity  u 


The  HEG  operating  conditions  I  -  IV  are  the  original  high  enthalpy  conditions  covering  a  total  specific 
enthalpy  range  from  12-23  MJ/kg.  A  detailed  discussion  of  these  conditions  is  given  in  Hannemann  et  al. 
[28].  They  have  e.g.  been  used  for  the  investigation  of  several  re-entry  configurations  including  ARD, 
X-38  [32]  or  Pre-X  /  IXV. 

Over  the  last  years  the  HEG  operating  range  was  subsequently  extended.  In  this  framework  the  main 
emphasis  was  to  generate  test  section  conditions  which  allow  investigating  the  flow  past  hypersonic 
configurations  at  low  altitude  Mach  6  up  to  Mach  10  flight  conditions  in  approximately  33  km  altitude. 
These  low  enthalpy  conditions  cover  the  range  of  total  specific  enthalpies  from  1.5  -6  MJ/kg.  For  1:1 
scale  wind  tunnel  models,  conditions  XIII  and  XIV  duplicate  M  =  7.4  flight  conditions  in  28  km  and 
33  km,  respectively.  They  were  used  for  the  ground  based  testing  of  the  HyShot  II  [53]  and  IV  supersonic 
combustion  flight  experiment  configurations.  Condition  XXI  is  related  to  M  =  6  flight  conditions  in  15  km 
flight  altitude  and  condition  XXII  represents  the  HEG  operating  condition  with  the  highest  unit  Reynolds 
number.  Condition  XXI  was  e.g.  used  in  the  framework  of  the  SHEFEX  I  post  flight  analysis  and  the 
investigation  of  the  intake  of  the  LAPCAT  (EC  project  “Long  Term  Advanced  Propulsion  Concepts  and 
Technologies”)  M  =  8  aircraft.  Condition  XXXI  focuses  on  M  =  10  flight  conditions  in  33  km  altitude.  It 
will  be  used  in  the  framework  of  the  DLR  SHEFEX  II  flight  experiment  and  was  applied  for  the  ground 
based  testing  of  a  1:1  scale  scramjet  flight  experiment  configuration. 

With  reducing  flight  velocity  or  i.e.  total  specific  enthalpy,  the  binary  scaling  parameter  becomes  less 
important.  However,  for  1:1  scale  models  of  hypersonic  flight  configurations,  the  pL  versus  u  diagram 
gives  an  indication  of  the  dynamic  pressure  range  which  can  be  reproduced  in  HEG.  In  particular  for 
airbreathing  propulsion  driven  hypersonic  vehicles  this  quantity  is  of  interest.  Typically  they  fly  at 
dynamic  pressures  between  50  and  100  kPa  and  in  order  to  model  the  combustion  process  in  e.g.  a 
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scramjet  engine  correctly,  the  dynamic  pressure  is  one  of  the  quantities  which  must  be  duplicated  [42].  In 
Figure  7,  the  low  enthalpy  HEG  operating  conditions  are  also  given  in  terms  of  Mach  and  Reynolds 
number.  As  reference  length,  a  fixed  wind  tunnel  model  length  of  0.3  m  and  the  lengths  of  the  considered 
flight  configuration  were  used. 


Figure  7:  HEG  operating  conditions  in  terms  of  Mach  and  Reynolds  number 


Details  of  the  HEG  operating  conditions  discussed  above  are  provided  in  Table  1.  The  test  time  for  the 
high  enthalpy  conditions  I  -  IV  is  up  to  1  ms.  For  the  low  enthalpy  conditions,  the  test  time  ranges  from 
3-6  ms. 


Nozzle 

2 

3 

4 

5 

Condition 

I 

II 

III 

IV 

XIII 

XIV 

XXI 

XXII 

XXXI 

P0  [MPa] 

35 

85 

44 

90 

17 

8 

37 

54 

70.0 

T„[K] 

9100 

9900 

7000 

8100 

2740 

2810 

1640 

1200 

4400 

h0  [MJ/kg] 

22 

23 

12 

15 

3.3 

3.4 

1.5 

1.3 

6.0 

Mx  [-] 

8.2 

7.8 

8.1 

7.9 

7.4 

7.4 

6.0 

6.1 

10.3 

Rem  [1/m- 106] 

0.20 

0.42 

0.39 

0.67 

3.70 

1.60 

45.0 

100.0 

2.8 

Px  [Pa] 

660 

1700 

790 

1680 

1990 

880 

20100 

29400 

930 

TJK] 

1140 

1450 

800 

1060 

266 

277 

221 

152 

253 

Pjg/m3] 

1.7 

3.5 

3.3 

5.3 

25.9 

11.0 

327.0 

682.0 

12.6 

uy  [m/s] 

5900 

6200 

4700 

5200 

2410 

2450 

1750 

1510 

3270 

Table  1 :  Summary  of  HEG  nozzle  reservoir  and  test  section  flow  conditions 
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In  order  to  realise  the  operating  conditions  of  Table  1,  a  series  of  different  Laval  nozzles  had  to  be 
designed,  constructed  and  implemented  in  HEG.  The  nozzles  used  to  generate  the  corresponding  test 
conditions  are  given  in  Table  1.  Details  of  the  four  operational  HEG  nozzles  are  provided  in  the  following 
section.  Further,  different  pistons  are  utilised  on  HEG  for  generating  different  operating  conditions.  In 
order  to  allow  a  large  flexibility  in  tuning  new  operating  conditions,  four  pistons  (without  brakes)  with 
different  weight  (275  kg,  481  kg,  700  kg  and  815  kg)  are  available.  One  additional  848  kg  piston  with 
brakes  is  utilised. 

2.2.2  Set  of  HEG  Nozzles 

The  set  of  HEG  nozzles  comprises  of  a  conical  nozzle  used  for  the  high  enthalpy  conditions  I-IV  and  three 
contoured  nozzles  for  the  low  enthalpy  conditions.  Their  nominal  design  Mach  number,  area  ratio  and 
length  are  given  in  Figure  7.  Please  note  that  for  nozzle  2  the  Mach  number  is  lower  than  the 
corresponding  flight  Mach  number  (see  also  Figure  6)  due  to  chemical  and  thermal  freezing  effects  during 
the  nozzle  expansion  [28].  However,  for  high  enthalpy  testing  the  Mach  number  is  of  less  importance  and 
the  flight  velocity  must  be  reproduced  correctly. 

Due  to  the  different  nozzle  length,  a  second  test  section  was  built  for  nozzle  3.  When  utilising  nozzle  4,  an 
additional  adapter  ring  is  used  between  the  second  and  the  main  test  section.  The  length  and  diameter  of 
the  original  test  section  is  1.6  m  and  1.2  m,  respectively.  The  second  test  section  is  a  copy  of  the  original 
test  section  with  slightly  reduced  dimensions  (1.0  m  length  and  1.0  m  diameter).  The  adapter  ring  has  a 
length  of  0.8  m  and  a  diameter  1.0  m.  The  nozzle  -  test  section  assembly  using  the  four  HEG  nozzles  is 
shown  in  Figure  9.  Depending  on  the  used  operating  condition  and  the  angle  of  attack,  model 
configurations  with  a  typical  length  between  0.4  m  and  1.0  m  and  a  width  of  up  to  0.4  m  can  be  mounted 
in  the  test  sections.  In  case  the  major  emphasis  of  the  tunnel  testing  is  on  the  investigation  of  internal  flow 
paths  (e.g.  scramjet  combustors),  models  of  up  to  2.0  m  length  can  be  used.  The  weight  of  the  models  is 
typically  less  than  200  kg. 


Nozzle  2  (condition  MV)  Nozzle  3  (condition  XIII,  XIV) 


Nozzle  4  (condition  XXI,  XXII) 


Nozzle  5  (condition  XXXI) 


Nozzle 

2 

3 

4 

5 

nominal  M 

8 

8 

6 

10 

d*[mm] 

22 

40 

50 

26 

ymmi 

880 

590 

430 

880 

^Exi/A 

1600 

218 

74 

1146 

L  [mm] 

3780 

2710 

2015 

3780 

Figure  8:  Geometry  of  operational  HEG  nozzles 
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Figure  9:  HEG  nozzle  -  test  section  assembly 


2.2.3  Recent  HEG  Infrastructure  Upgrades 

The  data  acquisition  system  of  HEG  consists  of  a  total  of  180  channels.  Two  different  types  of  transient 
recorders  are  currently  used.  The  original  KRENZ  System  from  Eckelmann  Industrie  Automation  with  60 
channels  is  characterised  by  1  MHz  sampling  rate  and  0.1  Mega  samples  per  channel  at  14  Bit  resolution. 
The  new  SATURN  System  from  AMO  GmbH  provides  120  channels.  Its  sampling  rate  is  50  MHz  and  10 
Mega  samples  per  channel  can  be  recorded  at  16  Bit  resolution. 

A  gaseous  hydrogen  injection  system  was  installed  at  HEG  in  order  to  allow  the  delivery  of  hydrogen  fuel 
to  wind  tunnel  models  for  the  investigation  of  scramjet  combustion.  The  fuel  system  consists  of  a  12  mm 
diameter  and  38.4  m  long  Ludwieg  tube,  and  a  fast  acting  solenoid  valve.  The  maximum  filling  pressure 
of  the  Ludwieg  tube  is  15  MPa  and  it  can  deliver  a  pulse  of  fuel  with  constant  pressure  for  up  to  50  ms. 


Figure  10:  HEG  calibration  rake  -  design  drawing  (left)  and  rake  installed  in  the  test  section 

(right) 
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A  new  modular  cross  arm  calibration  rake  was  designed  and  manufactured  which  can  be  used  for  the 
detailed  calibration  of  the  test  section  free  stream  flows  generated  by  nozzles  2-5.  Pitot  pressure,  static 
pressure  and  stagnation  point  heat  transfer  gauges  can  be  mounted  on  this  rake  (see  Figure  10). 


3.0  MEASUREMENT  TECHNIQUES  FOR  SHORT  DURATION  HYPERSONIC 
GROUND  BASED  TESTING 

In  this  section  an  overview  of  measurement  techniques,  which  are  suitable  for  short  duration  hypersonic 
ground  based  testing,  is  given.  The  emphasis  is  put  on  those  techniques  which  are  applied  in  the  High 
Enthalpy  Shock  Tunnel  Gottingen,  HEG. 


3.1  Pressure  Measurements 

Surface  pressure  measurements  in  short  duration  facilities  require  different  approaches  compared  to 
continuously  running  facilities.  Because  the  measurement  time  is  short,  pressure  transducers  with  fast 
response  time  have  to  be  used.  Additionally,  the  susceptible  area  of  the  transducer  must  be  installed  close 
to  the  surface  to  minimise  the  filling  time  of  the  tubing  system  in  front  of  the  susceptible  area.  The  most 
commonly  used  pressure  gauges  in  short  duration  hypersonic  ground  based  testing  are  based  on  the 
piezoelectric  and  on  the  piezoresistive  effect.  The  piezoelectric  effect  leads  to  the  generation  of  a  voltage 
in  a  crystal  and  the  piezoresistive  effect  leads  to  an  increase  of  the  resistance  of  a  semiconductor  when  a 
pressure  load  is  applied.  Optical  fibre  pressure  transducers  based  on  different  principles  including  intensity 
modulation,  interferometry,  polarization  effects,  refractive  index  changes,  reflectometry  and  fibre  Bragg 
grating  are  currently  under  development.  A  review  of  micromachined  pressure  gauges  is  given  e.g.  in  [14]. 
The  development  of  an  optical  pressure  gauge  based  on  a  Fabry-Perot  interferometer  with  a  natural 
frequency  above  1  MHz,  enables  data  capture  over  a  bandwidth  exceeding  100  kHz  ([93]).  This  results  in 
response  times  which  are  suitable  for  short  duration  hypersonic  testing. 


implanted  piezo  resistors 

Si02 


glas  body 


pressure  applied 


Figure  11 :  Schematic  of  a  Kulite  piezoresistive  pressure  transducer 


Miniaturised  fast  piezoresistive  pressure  transducers  are  e.g.  manufactured  by  Kulite.  Depending  on  the 
pressure  range,  the  Kulite  XCEL-100  pressure  gauge  with  a  diameter  of  2.4  mm  has  a  natural  frequency 
between  240  kHz  and  1  MHz  and  a  bandwidth  of  20  kHz.  Figure  1 1  shows  the  schematic  of  a  Kulite 
pressure  transducer.  A  silicon  element  is  used  as  mechanical  diaphragm  and  the  sensing  element  itself  is 
an  integral  part  of  this  silicon  element.  The  piezoresistors  are  formed  within  the  silicon  diaphragm  by 
either  diffusion  or  implantation  of  atoms.  Two  of  these  resistors  are  positioned  on  the  silicon  diaphragm 
such  that  they  experience  a  compressive  strain,  and  two  are  positioned  such  that  they  experience  a  tensile 
strain.  They  are  connected  to  form  a  fully  active  Wheatstone  bridge.  These  transducers  may  be 
manufactured  very  small-sized:  diameters  of  the  housing  can  be  in  the  range  of  1  to  2  mm. 
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To  achieve  fast  response  times,  it  is  necessary  to  position  the  susceptible  area  of  the  transducer  as  close  as 
possible  to  the  model  surface.  As  shown  in  Figure  12,  small  holes  are  drilled  in  the  model  which  act  as 
pressure  tappings.  The  gauge  is  placed  directly  behind  these  tappings.  In  Figure  12,  the  schematic  of  two 
different  types  of  pressure  gauge  installations  as  used  in  HEG  is  shown. 


^flow 


Figure  12:  Installation  of  pressure  transducers  into  a  model  wall  (transducers  shown  in  blue) 


Typical  tapping  diameters  range  from  0.5  to  1.5  mm.  The  schematic  on  the  left  side  of  Figure  12  shows  an 
installation  suited  for  tappings  perpendicular  to  the  flow.  This  installation  is  not  suited  for  stagnation 
regions,  where  the  flow  may  directly  stagnate  on  the  pressure  transducer,  leading  to  a  high  heat  flux  on  the 
susceptible  area.  A  stagnation  point  installation  is  shown  on  the  right  side  of  Figure  12.  Here  the 
susceptible  area  of  the  transducer  is  protected  by  a  stagnation  plate.  Examples  of  pressure  measurements 
using  both  configurations  are  given  in  Figure  13.  One  pressure  transducer  was  mounted  in  the  side  wall  of 
a  shock  tube  and  two  transducers  were  installed  in  the  end  wall,  measuring  the  pressure  rise  caused  by  the 
reflection  of  the  incident  shock  wave.  It  is  clear  that  the  stagnation  plate  leads  to  an  increase  of  the  rise 
time.  For  the  applied  tapping  diameter  and  cavity  size,  the  resulting  increase  amounts  to  about  0.15  ms. 


t3 


Figure  13:  Rise  time  investigation  for  different  pressure  gauge  installations  in  shock  tube  walls 


3.2  Surface  Heat  Flux  Measurements 

In  short  duration  high  enthalpy,  hypersonic  ground  based  test  facilities  the  most  common  way  to  measure 
wall  heat  flux  is  by  utilising  thin  film  gauges  or  thermocouples.  Due  to  the  harsh  environment  in  these 
facilities  with  high  wall  shear  stresses  during  the  starting  process  of  the  flow,  the  most  robust  and  reliable 
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technique  to  measure  heat  flux  is  the  application  of  surface  mounted  coaxial  thermocouples.  The  principle 
of  operation  of  thermocouples  is  the  Peltier-Seebeck  or  thermoelectric  effect  ([48]). 


Figure  14:  Schematic  of  MedTherm  type  E  thermocouple  (left)  and  photograph  (right) 


In  order  to  allow  a  detailed  instrumentation  of  wind  tunnel  models,  miniaturized  thermocouples  with 
diameters  as  small  as  0.4  mm  are  used.  Figure  14  shows  a  schematic  view  of  a  miniaturized  coaxial  ultra 
fast  thermocouple  manufactured  by  Medtherm  Corporation  ([56]).  Coaxial  thermocouples  are  built  with  a 
core  consisting  of  one  conductive  material  separated  by  insulation  from  the  outer  ring  made  of  a  different 
material.  Possible  material  combinations  and  their  official  ISA  (Instrument  Society  of  America)  code  are 
given  in  Table  2  together  with  their  application  range.  The  electrical  connection  on  the  top  surface  is 
realized  either  by  coating  the  surface  with  one  of  the  materials  or  simply  by  sanding  of  the  transducer 
surface.  In  either  case  the  mass  of  the  electrical  connection  in  comparison  to  the  mass  of  the  transducer 
body  is  small  leading  to  response  times  of  thermocouples  in  the  order  of  a  few  ps.  Additionally,  the  heat 
flux  into  the  transducers  is  dominated  by  the  thermocouple  body  and  the  junction  itself  is  negligible  due  to 
its  low  mass.  The  sensitivity  of  thermocouples  range  from  1  pV/K  to  70  pV/K.  The  thermocouple  inherits 
the  advantage,  that  the  shape  of  the  transducer  can  be  adapted  to  any  wall  curvature  by  sanding.  Further,  in 
case  of  transducer  failure,  it  can  also  be  repaired  in  situ  by  sanding  the  surface. 


ISA  code 

Material  A  (+) 

Material  B(-) 

Application  Range 

B 

Ratinum  30%  Rhodium 

Ratinum  6%  Rhodium 

1 640K  to  1970K 

C 

95%  W5Fte  Tungsten  5%  Rhenium 

W26Fte  Tungsten  26%  Rhenium 

1 920K  to  2585K 

E 

Chromel 

Con  st  ant  an 

365K  to  1170 

J 

Iron 

Con  st  ant  an 

365K  to  1030K 

K 

Chromel 

Alumel 

365K  to  1530K 

N 

Nicrosil 

Nisil 

920K  to  1530K 

R 

Ratinum  13%  Rhodium 

Ratinum 

1 1 40K  to  1720K 

S 

Ratinum  10%  Rhodium 

Ratinum 

1 250K  to  1720K 

T 

Copper 

Const  ant  an 

70K  to  620K 

Chromef 

Alumef 

Con st  ant  an 

Nisil+ 

NicrosiT 

90%  Nickel  10%  Chromium 

95%  Nickel  2%  Manganese  2%  Alumimnium  1%  silicon 

55%  Copper  45%  Nickel 

95.6%  Nickel  4.4%  Slicon 

84.1%  Nickel  1.3%  Slicon  14.6%  Chromium 

Trademarks  of  Hoskins  Manufacturing  Company,  +Trademarksof  Harrison  Alloys,  Inc. 

Table  2:  Possible  material  combinations  and  application  ranges  of  thermocouples 

Alternatively  to  thermocouples,  in  areas  with  less  harsh  flow  environments  such  as  wake  flows  or  surfaces 
at  small  angles  of  attack,  thin  film  gauges  can  be  used  to  measure  the  temperature  history  in  order  to 
evaluate  the  heat  flux  (e.g.  [62]  or  [60]).  Thin  film  gauges  consist  of  a  nonconductive  substrate  on  which  a 
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metallic  film  with  a  typical  thickness  of  less  than  1  pm  is  metallised.  The  working  principle  of  these 
gauges  is  based  on  the  dependence  of  the  electrical  resistance  of  the  metallic  film  on  temperature.  For  the 
metallic  film  noble  metals  are  used  because  they  exhibit  the  largest  sensitivity.  Again,  the  mass  of  the 
susceptible  area  is  very  small  compared  to  the  body  in  which  the  heat  is  conducted.  Therefore,  similar 
response  times  as  for  thermocouples  are  obtained. 


Figure  15:  Schematic  of  ISL  thin  film  gauge  (left)  and  photograph  (right);  (diagram  and 

photograph  courtesy  ISL) 

As  an  example,  in  Figure  15,  a  thin  film  gauge  as  manufactured  by  the  French-German  Research  Institute 
of  Saint-Louis  (ISL)  is  shown.  The  measuring  point  is  a  thin  platinum  film,  which  is  metallised  on  a  glass 
body.  The  electrical  connection  to  the  measurement  wires  is  realized  by  a  silver  coating  of  the  upper  and 
lower  third  of  the  glass  body.  The  typical  diameter  of  these  gauges  is  in  the  order  of  millimetres.  The 
temperature  sensitivity  and  signal  to  noise  ratio  of  thin  film  gauges  are  superior  to  thermocouples. 
However,  disadvantages  are  the  missing  robustness  in  harsh  flow  environment  as  usually  found  in  short 
duration  facilities,  the  large  effort  of  adapting  the  susceptible  area  to  curved  surfaces  and  the  missing 
possibility  to  repair  the  gauge  in  situ. 

Two  assumptions  are  used  to  evaluate  the  heat  flux  from  the  measured  temperature  traces  using  coaxial 
thermocouples  or  thin  film  gauges.  The  first  assumption  is  that  the  measurement  is  dominated  by  one 
dimensional  heat  conduction  in  the  sensor  body,  and  the  second  is  that  the  gauge  itself  acts  as  a  semi 
infinite  body  ([77]).  To  ensure  that  the  first  assumption  is  valid,  only  the  transducer  face  should  be 
exposed  to  the  flow,  and  the  combination  of  wind  tunnel  model  wall  material  and  transducer  material  has 

to  be  chosen  in  an  appropriate  way.  The  validity  of  the  semi  infinite  body  assumption  requires  v > \Ja-t 
([15]),  where  v  is  the  length  of  the  transducer,  a  the  diffusivity  and  t  the  time.  The  diffusivity,  a ,  is 
defined  as  the  ratio  of  the  heat  conductivity,  k  ,  and  the  product  of  density,  p ,  and  specific  heat  capacity, 
c  ,  a  =  k /  p-c .  For  thermocouples  of  ISA  code  E  with  a  length  of  10  mm  the  time  before  the  heat  reaches 
the  end  of  the  transducer  is  approximately  3  s.  This  period  of  time  is  orders  of  magnitude  higher  than  the 
typical  test  time  in  short  duration  facilities  of  a  few  milliseconds.  If  the  two  assumptions  mentioned  above 
can  be  regarded  as  valid,  the  heat  conduction  problem  in  the  transducer  body  is  described  by  the  following 
differential  equation: 


d2T  _  1  8T 
8x2  a  8t 

Here,  T  is  the  temperature,  x  the  coordinate  normal  to  the  wind  tunnel  model  wall  and  t  the  time.  With  the 
boundary  conditions  that  the  heat  flux  into  the  model  wall,  qw  (x  =  0)  =  —ic(dT  /  dx)  and  T(x  =  oo)  =  0 , 

the  following  equation  can  be  derived  for  the  evaluation  of  the  heat  flux  as  shown  by  Schultz  &  Jones 
[77]: 
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For  the  evaluation  of  the  wall  heat  flux  by  numerical  integration  of  digitally  stored  temperature  data,  the 
above  relationship  is  replaced  by  a  discretised  form  such  as 


4wO„)  =  2 


PCK , 


TC 


7=1 


+ 

1 

IL^ 

It  —t  , 

V  "  J  \ 

1 n  j~1 

(3) 


proposed  by  Cook  and  Felderman  [13]  or  an  alternative  scheme  described  by  e.g.  Kendall  et  al.  [44]. 


t  [ms]  after  shock  reflection  t  [ms]  after  shock  reflection 


Figure  16:  Heat  flux  evaluation  of  a  temperature  recording  in  a  shock  tunnel  (red)  compared  with 
the  temperature  evolution  based  on  a  stepwise  heat  transfer  load  (black) 

In  Figure  16  an  example  obtained  with  the  data  evaluation  method  of  Cook  and  Felderman  [13]  is 
presented.  On  the  left  side  of  Figure  16,  the  temperature  development  measured  by  a  type  E  thermocouple 
installed  in  the  stagnation  point  of  a  sphere  positioned  in  a  hypersonic  flow  with  total  specific  enthalpy  of 
h0  =  15  MJ/kg  is  indicated  by  the  red  line.  On  the  right  hand  side  of  Figure  16,  the  resulting  heat  flux  over 
time  (red  line)  is  shown.  If  the  heat  flux  into  the  wall  is  given  by  an  ideal  step  function,  the  exact  solution 
for  the  temperature  rise  at  x  =  0  resulting  from  equation  (1)  is: 


rw(*  =  0,f) 


I  t  = 

4tT  ]/  p'C'K 


(4) 


Using  a  constant  heat  flux  into  the  transducer  body  of  qw  =10  MW/m2  as  approximation  of  the  measured 
heat  flux  evolution  (the  black  line  in  the  right  part  of  Figure  16),  the  parabolic  temperature  history  profile 
shown  as  black  line  in  the  left  part  of  Figure  16  is  obtained. 


To  ensure  that  the  measurement  is  not  affected  by  conduction  in  wall  tangential  direction,  the 
thermocouple  type  should  be  chosen  in  such  a  way,  that  the  material  properties  of  the  model  wall  match 
those  of  the  thermocouple  at  the  measurement  position.  The  influence  of  different  thermocouple  /  model 
wall  material  combinations  on  the  measured  temperature  distribution  is  shown  in  Figure  17  and  Figure  18. 
These  results  were  obtained  by  numerical  investigation  of  the  unsteady  heat  conduction  process  in  a 
type  E  thermocouple  with  1.5  mm  diameter  and  a  length  of  15  mm,  installed  in  model  walls  of  different 
materials,  using  Ansys  (CAD-FEM  GmbH).  The  thermocouple  is  modelled  as  a  solid  cylinder  consisting 
of  chromel.  At  the  top  surface  a  constant  heat  flux  of  qw  =  10  MW/m2  is  applied.  The  temperature 
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evolution  at  the  contact  point  between  the  thermocouple  and  the  model  wall  is  given  in  Figure  17  for  three 
different  model  wall  materials.  Dependent  on  the  material,  large  differences  of  the  temperature  evolution 
are  obtained. 


qw=10MW/m2 


1400  [Kl 

300  340  S00  1000 


Figure  17:  Unsteady  finite  element  analysis  of  the  heat  conduction  process  in  a  type  E 
thermocouple  and  the  adjacent  model  wall  for  a  wall  heat  flux  of  10  MW/m2;  axisymmetric 
temperature  distribution  at  t=4  ms  for  a  carbon  fibre  model  wall  (left)  and  time  evolution  of  the 
temperature  at  the  interface  (A)  between  the  thermocouple  and  the  model  wall  for  different  wall 

materials  (right) 


y  [mm] 


y  [mm] 


Figure  18:  Difference  of  radial  computed  temperature  distribution  at  the  model  wall  and 
thermocouple  surface  obtained  for  different  thermocouple  /  wall  material  combinations 
compared  to  the  combination  chromel  /  chromel  for  which  no  temperature  difference  between 

thermocouple  and  wall  exists 


The  difference  of  the  radial  temperature  distribution  at  the  model  wall  and  thermocouple  surface  obtained 
for  the  thermocouple  /  wall  material  combination  chromel  /  chromel,  i.e.  a  combination  for  which  no 
temperature  difference  between  the  thermocouple  and  the  model  wall  exists,  and  the  material 
combinations  chromel  /  steel,  chromel  /  aluminium  and  chromel  /  carbon  fibre  are  plotted  in  Figure  18.  On 
the  left  side  of  Figure  18,  a  zoom  of  the  core  region  of  the  thermocouple  is  shown.  The  right  side  of  Figure 
18  shows  the  distribution  up  to  a  radius  of  1  mm.  In  case  of  a  coated  thermocouple  with  dimensions  as 
shown  in  Figure  14,  the  temperature  measurement  takes  place  at  y  =  0  -  0.0625  mm,  the  location  of  the 
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chromel  /  constantan  junction.  For  a  thermocouple  which  is  adapted  to  a  three  dimensional  surface  or 
repaired  in  situ  by  sanding,  the  chromel  /  constantan  junction  is  not  precisely  defined  and  the  complete 
thermocouple  surface  may  act  as  the  measurement  location.  As  can  be  seen  from  Figure  18,  depending  on 
the  thermocouple  /  wall  material  combination  the  assumption  of  one  dimensional  heat  conduction  can  be 
violated. 

From  equation  (2)  it  becomes  evident,  that  the  error  in  the  estimation  of  the  heat  flux  is  directly 
proportional  to  the  error  resulting  from  the  determination  of  yj pcK  .  To  accurately  measure  this  value  two 

procedures  may  be  applied.  The  contact  procedure  is  used  for  thermocouples  and  thin  film  gauges  and  is 
based  on  the  fact  that  under  the  assumption  of  one  dimensional  heat  conduction,  the  contact  temperature  of 

two  media  is  only  dependent  on  their  temperature  before  contact  and  on  the  thermal  properties  ( yj pcK  )  of 

both  media.  Knowing  the  thermal  properties  of  one  media  allows  to  determine  the  gauge  properties.  A 
possible  setup  is  a  mechanically  driven  device  that  generates  an  instantaneous  contact  between  the  gauge 
and  a  liquid  surface.  Such  a  method  is  discussed  by  lessen  and  Gronig  [41].  The  electric  discharge 
calibration  can  only  be  applied  for  thin  film  gauges.  It  utilizes  the  effect  that  with  an  electric  discharge,  the 
thin  film  itself  can  be  heated  with  a  known  and  constant  heat  flux.  Following  equation  (3),  the 
measurement  of  the  temperature  of  the  thin  film  gauge  results  in  the  determination  of  pcK  .This  procedure 
is  described  by  e.g.  Skinner  [81]. 

3.3  Force  and  Moment  Measurements 

The  ability  to  perform  integrated  force  and  moment  measurements  in  ground  based  testing  facilities  is  an 
important  part  of  the  design  and  development  of  hypersonic  vehicles.  For  a  force  balance  to  operate  in 
short  duration  ground  based  test  facilities  with  test  times  of  the  order  of  milliseconds  or  less,  it  would 
become  necessary  for  the  balance  to  have  an  extremely  short  response  time.  However,  given  the  short  test 
times  and  the  response  times  of  force  measurement  techniques  for  convential  wind  tunnels  (see  e.g.  Ewald 
[17]),  static  equilibrium  between  the  model  and  support  structure  is  rarely  established,  i.e.  it  may  only  be 
obtained  using  unrealistically  small  models.  Therefore,  it  becomes  necessary  to  use  specialised 
measurement  methods  that  account  for  the  dynamic  response  of  the  system.  An  overview  of  the  basic 
force  measurement  techniques  that  have  been  applied  in  short  duration  flows  is  given  e.g.  by  Robinson 
[72]. 

In  HEG  a  stress  wave  force  balance  was  implemented  [73].  This  technique  relies  on  the  ability  to  measure 
the  dynamic  response  of  the  model  and  supporting  structure.  Consequently,  any  effects  such  as  model 
flexibility  and  mass  distribution  are  accounted  for  (see  e.g.  Robinson  [72]).  Stress  wave  force  balances 
require  a  dynamic  calibration  of  the  complete  model,  balance  and  support  structure  assembly. 

Upon  flow  arrival,  stress  waves  propagate  through  the  model  at  the  speed  of  sound  of  the  material  and 
subsequently  enter  a  stress  bar  that  is  instrumented  with  a  strain  gauge  to  record  the  time  history  of  strain. 
If  the  system  is  linear,  the  resulting  strain  signal,  y(t),  can  be  related  to  the  applied  aerodynamic  load,  u(t), 
via  an  impulse  response  function,  g(t),  as  described  by  the  convolution  integral, 


y(t)  =  jg(t~T)u(T)dT.  (5) 

0 

The  aerodynamic  force  in  an  experiment  can  be  determined  by  the  deconvolution  of  the  strain  signal  with 
the  impulse  response  function.  The  impulse  response  function  is  determined  either  through  experimental 
calibration  or  through  finite  element  analysis.  However,  in  order  to  reduce  errors  due  to  modelling 
approximations  it  is  usually  preferable  to  determine  the  impulse  response  function  experimentally. 
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Sanderson  and  Simmons  [75]  first  demonstrated  the  technique  by  measuring  the  drag  force  on  a  15° 
conical  model,  200  mm  in  length  in  a  reflected  shock  tunnel.  The  model  was  made  from  aluminium  and 
attached  to  a  2  m  long  hollow  brass  bar.  This  stress  bar  or  sting,  was  instrumented  with  strain  gauges  to 
record  the  axial  strain  time  history.  The  test  time  was  approximately  1  ms  and  the  aerodynamic  load  was 
subsequently  found  through  a  numerical  deconvolution  process.  The  measured  drag  was  found  to  be  in 
good  agreement  with  theoretical  calculations. 

The  technique  was  then  extended  to  measure  multiple  components  of  force  by  Mee  et  al.  [57].  This  was 
done  on  a  15°  cone  with  an  internal  balance  arrangement.  The  lift,  drag  and  pitching  moments  were 
measured  by  short,  stiff  stress  bars  instrumented  with  strain  gauges.  The  technique  has  subsequently  been 
applied  to  a  wide  variety  of  models,  including  models  with  simultaneous  measurements  of  pressure  and 
heat  flux.  The  major  disadvantage  with  this  technique  is  that  aerodynamic  shielding  and  vibration  isolation 
of  the  support  structure  is  required  in  order  to  separate  the  stress  waves  generated  in  the  model  from  that 
of  the  test  section  environment.  Calibration  of  the  model  is  usually  performed  using  an  instrumented 
impact  hammer  or  via  a  cut-weight  technique.  In  order  to  maximise  the  performance  of  stress  wave  force 
balances,  an  individual  design  for  each  configuration  to  be  tested  is  recommended. 

Stress  wave  force  balances  have  been  developed  with  external  as  well  as  internal  stress  bar  arrangements. 
The  internal  three  component  force  balance  applied  by  Robinson  and  Hannemann  [73]  in  HEG,  is  shown 
in  Figure  19.  Four  short  stiff  stress  bars  are  mounted  on  a  sting  and  each  bar  is  instrumented  with 
semiconductor  strain  gauges  to  measure  the  time  history  of  strain.  The  balance  is  able  to  measure  forces 
(approximately  50  to  5000  N)  within  1  -  5  ms  on  instrumented  models  at  angles  of  attack  from  -40  -  20°. 
The  accuracy  of  the  force  balance  is  estimated  at  approximately  ±5%  for  the  axial  component  and  ±4%  for 
the  normal  and  pitching  moment  components. 

Additional  force  measurement  techniques  based  on  external  stress  wave  force  balances,  accelerometer 
based  and  free  flight  based  force  measurement  techniques  are  currently  under  development  in  HEG. 


sting 


Figure  19:  Finite  element  discretisation  of  the  internal  stress  wave  force  balance  -  cone  model 

set  up  used  in  the  HEG  shock  tunnel 

3.4  Optical  Measurement  Techniques 
3.4.1  Phase  Step  Holographic  Interferometry 

Interferometry  may  be  used  as  a  technique  to  measure  the  variation  of  the  refractive  index  in  a  gaseous 
flow  in  the  test  section  of  a  short  duration  facility.  This  information  can  be  used  to  evaluate  the  density 


8-18 


RTO-EN-AVT-1 86 


NATO 

OTAN 


A  Closely  Coupled  Experimental  &  Numerical  Approach  for  Hypersonic  &  High 
Enthalpy  Flow  Investigations  Utilising  the  HEG  Shock  Tunnel  &  the  DLR  TAU  Code 


distribution  of  the  investigated  flow  field.  In  particular,  holographic  interferometry  does  not  require 
machining  or  manufacturing  of  test  section  windows,  mirrors  or  lenses  with  high  precision,  because 
parasitic  effects  caused  by  imperfections  of  these  components  cancel  out  when  applying  the  holographic 
two  step  procedure.  Therefore,  related  to  the  application  in  short  duration  ground  based  test  facilities,  this 
technique  has  replaced  the  classical  and  labour  intensive  Mach-Zehnder  interferometry.  In  the  subsequent 
paragraphs,  a  brief  introduction  of  the  phase  step  holographic  interferometry  technique  is  given  and  the  set 
up  used  at  HEG  is  described. 

The  absolute  speed,  c0 ,  at  which  light  propagates  in  vacuum  is  a  constant.  In  any  kind  of  gaseous  media 
the  speed  of  light,  c  ,  will  be  lower.  The  ratio  of  the  two  speeds  defines  the  refractive  index: 


n(p)  =  ^  =  \  +  Kxp.  (6) 

c 

The  Gladstone-Dale  relation  describes  that  in  a  gaseous  media,  consisting  of  one  species,  the  refractive 
index  depends  on  the  density,  p ,  and  the  Gladstone-Dale  constant  Kx  (see  e.g.  Merzkirch  [59]).  The 
latter  is  weakly  dependent  on  the  wavelength  and  is  specific  for  each  gas.  For  gas  mixtures  the  refractive 
index  is  given  by: 


n(p)=l  +  piX£,  (7) 

i= 1 

where  K?:  are  the  Gladstone-Dale  coefficients  for  the  single  gas  species,  are  the  species  mass  fractions, 

S  is  the  number  of  species  and  X  is  the  wavelength  of  the  laser  light  source.  The  definition  of  a  linearly 
composed  Gladstone-Dale  constant  applies  not  only  to  air  and  other  neutral  gas  mixtures  but  also,  in  high 
temperature  gas  dynamics,  to  chemically  homogenous  gases  where  the  molecules  are  either  in  different 
excited  states,  dissociated,  or  even  ionized. 


Figure  20:  Schematic  of  interference  experiment 

The  basic  principle  of  interferometry  is  sketched  in  Figure  20.  Two  rays  of  coherent  light  interfere  in  point 
P.  Each  ray  passes  a  zone  with  different  refractive  index,  which  leads  to  a  time  shift,  At ,  due  to  the 
different  propagation  speeds  of  the  light  in  the  two  zones: 

At  - - =  —  (n2  -  nx ) .  (8) 

C2  Cl  Co 
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If  the  resulting  difference  in  the  optical  path  length  equals  the  wavelength,  the  phase  shift  (j)  between  both 
rays  is  equal  to  2 n  .  Therefore,  the  following  relationship  can  be  derived: 

-^-  =  j(n2-nl)  =  ^(p2-pl).  (9) 

The  light  intensity  measured  in  point  P  is  proportional  to  the  phase  shift  between  the  two  rays. 
Consequently,  the  measurement  of  the  intensity  or  phase  shift  in  point  P  is  directly  related  to  the  density 
difference  in  regions  1  and  2  in  Figure  20. 


Figure  21 :  Principle  of  Mach-Zehnder  (left)  and  holographic  interferometry  (middle,  right) 


Figure  22:  Schematic  setup  of  the  HEG  holographic  interferometry  system 

Using  a  Mach-Zehnder  interferometer,  which  is  shown  schematically  on  the  left  side  of  Figure  21,  the 
density  in  the  test  section  can  be  evaluated  if  the  density  distribution  is  known  at  a  reference  point.  The 
schematic  emphasises  that  any  imperfection  which  leads  to  a  modification  of  the  light  path,  disturbs  the 
interference  measurement  in  plane  F.  To  avoid  this,  the  measurement  can  be  performed  in  two  steps  (see 
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middle  sketch  of  Figure  21).  In  the  first  step,  the  light  beam  (green)  passing  through  the  evacuated  test 
section  is  recorded  and  in  the  second  step  the  light  beam  (red)  passing  through  the  test  section  with  flow  is 
recorded.  Subsequently,  both  beams  are  reconstructed  resulting  in  an  interference  pattern  in  the 
measurement  plane,  P  (see  right  sketch  of  Figure  21).  The  advantage  of  this  technique  is  that  any 
imperfection  of  the  optical  setup  modifies  the  path  of  both  light  beams.  Therefore,  they  cancel  out  and  do 
not  influence  the  reconstructed  interference  measurement.  To  store  and  reconstruct  both  beams,  a 
holographic  storage  technique  has  to  be  used.  Therefore,  this  two  step  measurement  technique  is  called 
holographic  interferometry.  A  more  detailed  discussion  of  this  technique  can  be  found  e.g.  in  [89]. 

A  schematic  of  the  holographic  interferometry  setup  used  at  HEG  is  shown  in  Figure  22.  The  green  line 
shows  the  light  path  of  the  reference  beam  and  the  red  line  shows  the  light  path  of  the  object  beam.  The 
object  light  beam  passes  through  the  test  section  and  is  brought  to  interference  with  the  reference  beam  on 
the  holographic  plate.  Here  the  interference  pattern  between  both  beams  is  recorded.  To  achieve 
interference  between  the  object  and  the  reference  beam,  a  light  source  with  sufficient  coherence  length  is 
needed.  A  seeded  Innolas  Nd:YAG  Laser  (Model  Spitlight  300),  emitting  light  at  532  nm  is  used  because 
it  has  a  coherence  length  larger  than  1  m.  The  optical  path  length  of  one  of  the  beams  is  approximately 
15  m  and  the  successful  set  up  of  the  optical  system  is  facilitated  by  aligning  the  path  length  difference 
between  the  two  beams  to  the  coherence  length  of  the  laser. 

The  first  exposure  of  the  holographic  plate  with  one  reference  beam  is  performed  prior  to  a  run  in  HEG 
and  a  second  exposure  with  the  second  reference  beam  is  obtained  during  the  test  time.  After  the  chemical 
treatment  of  the  holographic  plate,  two  reconstruction  waves  are  created  in  a  separate  reconstruction  unit. 
These  reference  waves  used  in  the  reconstruction  unit  are  identical  to  the  reference  waves  used  for  both 
exposures. 

3.4.2  High  Speed  Flow  Visualization 

To  visualise  the  flow  during  experiments  in  HEG,  a  high  speed  flow  visualisation  system  is  utilised.  It 
uses  a  Z-path  layout  for  the  object  light  path  with  spherical  mirrors  (Halle  SDH4300)  having  a  diameter  of 
300  mm  and  a  focal  length  of  1500  mm.  A  schematic  of  the  optical  setup  is  shown  in  Figure  23.  Two 
different  light  sources  are  used.  The  diode  pulsed  Nd:YAG  laser  from  Lightwave  Electronics  (Model  612) 
emits  light  at  532  nm  and  has  the  capability  to  operate  at  pulse  rates  up  to  50  kHz.  Based  on  the 
specifications  of  the  manufacturer,  the  HSPS  (High-Speed  Photo-Systeme)  NANOLITE  KL-K  sparcflash 
lamp  allows  repetition  rates  up  to  20  kHz  with  fixed  pulse  duration  of  8  ns.  If  the  distance  of  the  electrodes 
is  reduced  and  the  spare  gap  is  operated  in  an  argon  atmosphere  instead  of  air,  the  repetition  rate  can  be 
increased  up  to  32  kHz.  The  light  is  widened  by  a  telescope  system  (LI  and  L2)  and  subsequently 
expanded  by  a  lens  (L3)  onto  the  spherical  main  mirror  (HI)  which  generates  a  parallel  light-bundle  of 
diameter  300  mm  through  the  test  section.  On  the  opposite  side  of  the  test  section  the  beam  is  collimated 
by  the  second  spherical  main  mirror  (H2)  onto  the  plane  mirror  S3.  With  lens  L4,  the  test  section  image  is 
focused  onto  the  film  plane  (A).  In  order  to  record  Schlieren  images,  a  razor  blade  (R)  is  positioned  at  the 
focal  point  of  the  spherical  mirror  H2.  For  flow  visualisation  using  the  shadowgraph  technique,  the  razor 
blade  is  removed  and  the  image  plane  in  the  test  section  is  slightly  shifted  by  moving  lens  L4  thus 
generating  a  shadowgraph  effect.  Image  recording  is  performed  with  two  different  cameras.  If  a  Cordin 
rotating  drum  camera  (Model  318)  is  used,  the  images  are  recorded  on  black  and  white  film  (Kodak 
TMAX  100).  The  camera  consists  of  a  rotating  drum,  which  is  equipped  with  35  mm  film  stripes  of  1  m 
length.  Alternatively,  a  Shimadzu  HPV-1  digital  high  speed  camera  can  be  applied.  This  camera  is  able  to 
record  up  to  100  frames  at  a  maximum  imaging  rate  of  1  MHz.  Using  the  HPV-1  camera,  the  effective 
exposure  time  is  defined  by  the  pulse  length  of  the  laser  of  approximately  60  ns  for  a  pulse  rate  of  16  kHz 
and  119  ns  for  32  kHz.  The  shutter  of  the  drum  camera  remains  open  during  the  experiments.  Therefore, 
the  pulse  rate  of  the  laser  defines  the  framing  rate  and  its  narrow  pulse  width  the  exposure  time.  Typical 
laser  pulse  rates  of  9  kHz,  10  kHz  and  15  kHz  result  in  exposure  times  of  approximately  45  ns,  48  ns  and 
56  ns,  respectively.  Therefore,  for  both  cameras,  the  high  speed  flow  visualisation  system  relies  on  the 
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short  duration  of  the  laser  light  pulse  to  record  a  still  image  during  the  significantly  larger  exposure  times 
of  the  cameras. 


Figure  23:  Schematic  of  the  optical  setup  for  the  high  speed  flow  visualisation  in  HEG.  H: 
Parabolic  mirrors,  S:  Plane  mirrors,  L:  Lenses,  A:  Focal  plane,  R:  Razor  blade 


4.0  NUMERICAL  TOOL  -  DLR  TAU  CODE 

4.1  General  Overview 

At  DLR,  the  development  and  utilisation  of  CFD  schemes  which  are  able  to  model  hypersonic  and  high 
enthalpy  flows  was  started  some  20  years  ago.  In  a  first  step,  the  physico-chemical  models  necessary  to 
compute  high  enthalpy  flows  were  implemented  in  the  implicit,  upwind,  Total  Variation  Diminishing 
(TVD)  Thin-Layer  Navier-Stokes  code,  NSHYP  [7], [8], [9]).  NSHYP  was  extensively  used  for  the 
computation  of  steady  two-dimensional  and  axisymmetric  high  enthalpy  flow  fields  (see  e.g.  [10], 
[28], [29]).  In  order  to  achieve  efficient  steady  state  high  enthalpy  flow  field  computations  past  complex 
three-dimensional  configurations  the  DLR  CEVCATS-N  code  was  developed  ([68], [11]).  The  physico¬ 
chemical  models  which  were  tested  in  NSHYP  were  included  into  the  DLR  CEVCATS  code,  a  block- 
structured,  three  dimensional  finite-volume  scheme  ([45] [46]).  It  uses  multigrid  strategies  with  residual 
averaging  and  local  time-stepping  to  accelerate  convergence  to  a  steady  state. 

Finally,  the  DLR  TAU-code  [21],  [79],  [71]  was  developed.  It  is  a  CFD  platform  for  the  computation  of 
viscous  and  inviscid  flows  which  is  suitable  for  the  investigation  of  steady  and  unsteady  flows  past 
complex  geometries.  It  covers  the  low  subsonic  up  to  the  hypersonic  flow  regime.  TAU  comprises  several 
means  for  grid  modification,  namely  the  adaptation  and  the  deformation  module.  Both  (block-)  structured 
and  hybrid  unstructured  grids  composed  of  hexahedrons,  prisms,  tetrahedrons  and  pyramids  can  be 
utilised.  The  first  two  element  types  of  unstructured  grids  are  usually  used  in  semi-structured  layers  above 
surfaces  allowing  for  a  better  resolution  of  boundary  layers.  Tetrahedrons  are  used  to  fill  the 
computational  domain  in  a  flexible  way,  allowing  for  local  refinement  without  hanging  nodes  while  the 
pyramids  are  needed  for  the  transition  between  elements  with  quadrilateral  faces  and  elements  with 
triangular  faces. 
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For  parallel  computations  the  grids  are  partitioned  into  a  given  number  of  domains  before  the  start  of  the 
computation.  Load  balancing  is  performed  on  edge-  and  point  weights  which  are  adjusted  to  the 
specifications  of  the  solver.  After  the  grid  partitioning,  all  other  modules  of  TAU  run  on  a  single  domain 
per  process.  Grid  re -partitioning  is  performed  either  if  the  grid  was  locally  (de-)refined  or  if  the  number  of 
domains  is  changed. 

The  standard  solver  module  uses  an  edge-based  dual-cell  approach  based  on  a  vertex-centered  scheme. 
The  pre-processing  module  computes  the  required  dual  grid  composed  of  general  control  volumes  from 
the  primary  elements.  They  are  stored  in  an  edge  based  data  structure,  which  makes  the  solver  independent 
of  the  element  types  of  the  primary  grid.  All  metrics  are  given  by  normal  vectors,  representing  size  and 
orientation  of  the  faces,  the  geometric  coordinates  of  the  grid  nodes  and  the  volumes  of  the  dual  grid  cells. 
The  connectivity  of  the  grid  is  given  by  linking  the  two  nodes  on  both  sides  of  each  face  of  the  dual  grid  to 
the  corresponding  edge  from  the  primary  grid  elements.  In  order  to  enable  the  use  of  a  multi-grid 
technique  an  agglomeration  approach  is  employed  to  obtain  coarse  grids  by  fusing  fine  grid  control 
volumes  together. 

In  the  solver  module,  inviscid  terms  are  computed  employing  either  a  second-order  central  scheme  or  a 
variety  of  upwind  schemes  using  linear  reconstruction  (of  the  left  and  right  states  of  a  dual  grid  face)  for 
second-order  spatial  accuracy.  Viscous  terms  are  generally  computed  with  a  second-order  central  scheme. 

Various  explicit  Runge-Kutta  schemes  and  an  implicit  approximate  factorization  scheme  (Lower-Upper 
Symmetric  Gauss-Seidel)  can  be  used  for  time  integration.  Additional  convergence  acceleration  is 
achieved  by  the  application  of  full  multi-grid  and  residual  smoothing  algorithms. 

For  time  accurate  computations  a  Jameson-type  dual  time  stepping  approach  [39]  is  employed.  Both,  grid 
deformation  as  well  as  bodies  in  arbitrary  motion  can  be  simulated  in  this  framework. 

The  RANS  turbulence  models  implemented  in  the  TAU  code  include  linear  as  well  as  non-linear  eddy 
viscosity  models  originating  from  both,  the  one-  as  well  as  the  two-equation  model  families.  The  standard 
turbulence  model  in  TAU  is  the  Spalart-Allmaras  model  yielding  highly  satisfactory  results  for  a  wide 
range  of  applications  while  being  numerically  robust.  Additionally,  a  number  of  two-equation  models 
based  on  a  k-oo  formulation  are  available.  Also  nonlinear  explicit  algebraic  Reynolds  stress  models 
(EARSM)  have  been  integrated.  The  implementation  and  validation  of  Reynolds  stress  models  is  ongoing 
work.  Further,  options  exist  to  perform  Detached  Eddy  Simulations  (DES)  based  on  the  Spalart-Allmaras 
[83]  or  the  Menter  SST  [58]  models  or  the  so-called  Extra-Large  Eddy  Simulation  (XLES).  In  order  to 
allow  the  modelling  of  transitional  flows,  the  turbulent  production  terms  are  suppressed  in  regions  which 
are  flagged  in  the  grid  as  being  of  laminar  flow  type. 

Detailed  descriptions  of  the  modelling  capabilities,  solution  algorithms  and  auxiliary  tools  of  the  TAU 
CFD  platform  are  given  in  [7 8], [79]  and  [20]. 

4.2  Modelling  of  compressible  and  reacting  high  enthalpy  flows 

The  DLR  Tau  code  includes  extensions  for  chemical  and  thermal  non-equilibrium  flows  in  high  enthalpy 
aerothermodynamics.  The  flow  is  considered  to  be  a  reacting  mixture  of  thermally  perfect  gases.  A 
dedicated  transport  equation  is  solved  for  each  individual  species.  The  chemical  source  term  in  this  set  of 
transport  equations  is  computed  from  the  law  of  mass  action  by  summation  over  all  participating  reactions. 
The  forward  reaction  rate  is  computed  from  the  modified  Arrhenius  law  and  the  backward  rate  is  obtained 
from  the  equilibrium  constant  which  is  computed  directly  from  the  partition  functions  of  the  participating 
species. 
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The  thermodynamic  properties  (energy,  entropy,  specific  heat)  are  calculated  from  the  partition  functions 
or  external  lookup  tables  for  each  individual  species  in  the  reacting  gas  mixture.  The  advantage  of  this 
approach  is  its  high  flexibility.  Extensions  such  as  multi  temperature  models  to  handle  thermal  non¬ 
equilibrium  effects  are  easily  possible.  Having  determined  the  mixture  composition  and  the 
thermodynamic  state  of  the  individual  species,  the  properties  of  the  reacting  gas  mixture  are  then 
computed  using  suitable  mixture  rules  such  as  proposed  by  Wilke  [95]  for  the  viscosity  and  by  Herning 
and  Zipperer  [35]  for  the  heat  conductivity. 

For  fully  catalytic  wall  boundaries,  a  Dirichlet  condition  for  the  species  mass  fractions  is  set  according  to 
the  local  equilibrium  composition.  Non-catalytic  walls  are  modeled  using  a  von  Neumann  boundary 
condition  imposing  vanishing  wall-normal  gradients  of  the  species  mass  fractions. 

The  species  diffusion  fluxes  are  modeled  using  Fick’s  law  applying  an  averaged  diffusion  coefficient  for 
all  species.  This  approximate  diffusion  coefficient  is  computed  using  the  mixture  viscosity  and  constant 
laminar  Schmidt  numbers,  Sc.  Turbulent  diffusion  is  modeled  in  an  analog  way  by  computing  a  turbulent 
diffusion  coefficient,  DT ,  from  the  eddy  viscosity,  and  the  turbulent  Schmidt  number,  ScT.  The  eddy 
viscosity  is  derived  from  the  applied  turbulence  model  (e.g.  computed  from  the  turbulent  kinetic  energy 
and  the  length  scale  when  applying  a  k-co  model). 

Thermal  non-equilibrium  flows  are  computed  by  solving  an  additional  transport  equation  for  the 
vibrational  energy  of  each  molecule  in  non-equilibrium.  The  relaxation  of  vibrational  energy  is  modeled 
according  to  the  Landau-Teller  [47]  approach  and  the  vibrational  relexation  times  are  obtained  from  the 
correlation  of  Millikan  and  White  [61]. 

An  assumed  Probability-Density-Function  (PDF)  model  as  described  by  Gerlinger  [24]  is  implemented  to 
model  the  influence  of  turbulent  fluctuations  on  the  species  source  terms  for  detailed  chemistry 
mechanisms.  Point  statistical  PDF  methods  in  conjunction  with  a  detailed  chemistry  scheme  have  a  very 
wide  range  of  application  and  have  been  successively  applied  to  e.g.  turbulent  supersonic  combustion 
phenomena  (Gerlinger  [24],  [22],  Gaffney  [19]).  There  are  no  limitations  concerning  their  applicability  for 
premixed  and  non-premixed  combustion  and  for  different  Damkohler  numbers.  The  averaged  turbulent 
chemical  source  term  is  computed  by  PDF-weighted  integration  over  the  parameter  space  (temperature, 
species  concentration).  Two  major  approaches  can  be  distinguished:  Assumed-PDF  methods  prescribe  the 
mathematical  shapes  and  the  PDF  function  is  usually  defined  by  its  1st  and  2nd  moments.  In  the  transport 
equation  PDF  approach  the  evolution  of  the  complete  PDF  function  is  computed  which  is  computationally 
much  more  expensive  but  offers  the  advantage  of  more  physical  PDF  shapes. 


4.2.1  Navier-Stokes  equations  of  a  mixture  of  reacting  gases 

The  Navier-Stokes  equations  for  a  reacting  mixture  of  compressible  ideal  gases  can  be  written  in  their 
integral  form  as: 


d_ 

dt 


J  UdV  +  \  FEu  ndS  =  \  FNS  ndS  +  J  QdV 

VS  S  V 


(10) 


The  vector  of  the  conservative  variables  in  the  case  of  thermal  equilibrium  (only  one  overall  energy 
equation  is  needed)  is: 


U=[ps,puT,pE^j 

The  matrix  of  the  inviscid  (Euler)  fluxes  is: 


(11) 
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And  the  matrix  of  the  viscous  (Navier-Stokes)  fluxes  reads: 
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Note  that  the  diffusion  flux,  psuds ,  of  species  S  is  modelled  using  Fick’s  law  and  an  averaged  diffusion 
coefficient  D  for  all  species: 


Pslls  =-pDV 
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This  approximate  diffusion  coefficient  is  computed  using  the  viscosity  and  a  constant  Schmidt  number,  Sc. 
The  viscous  stress  tensor,  P ,  is  modelled  using  the  Boussinesq  approximation: 


P  =  P 


VuT +(Vu)T  -|^(vrw)/ 


(15) 


The  source  vector  Q  includes  the  chemical  source  term  a>s : 


Q  = 


0)s 

0 

0 


(16) 


4.2.2  Computation  of  the  chemical  source  term 

In  addition  to  solving  the  set  of  transport  equations  that  describe  the  evolution  of  the  different  species  S  in 
the  reacting  flow  field,  the  computation  of  the  chemical  source  term  which  determines  the  rate  of 
production  or  destruction  of  a  species  due  to  chemical  reactions  is  a  key  part  of  the  reacting  flow  model. 
Chemical  sources  cos  are  produced  by  sets  of  chemical  reactions  involving  the  species  Xs  and  the 
stoichiometric  coefficients  as  and  J3S: 

I ysxs <  ^  >Zfisxs  (i7) 

s  s 

The  source  term  cos  can  be  computed  from  the  law  of  mass  action  by  summation  over  all  participating 
reactions: 
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The  forward  reaction  rate  is  obtained  from  the  modified  Arrhenius  law: 
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The  reaction  rate  of  the  backward  reaction  is  obtained  from  the  equilibrium  constant: 


k  =- 


k: 


(20) 


The  equilibrium  constant  can  be  directly  computed  from  the  partition  functions  of  the  participating 
species: 


Kq=U 


f  n  \&~a 

x! s 


(21) 


4.2.3  Thermodynamic  properties  and  partition  functions 

The  calculation  of  the  thermodynamic  properties  and  equilibrium  constants  for  the  individual  chemical 
reactions  is  based  on  either  partition  functions  or  external  lookup  tables.  The  partition  function  approach 
facilitates  the  application  of  multi  temperature  models  to  account  for  thermal  non-equilibrium. 

If  the  assumption  of  a  harmonic  oscillator  and  rigid  rotator  with  no  interdependencies  between  the  energy 
modes  is  made,  the  partition  function  of  a  species,  S,  can  be  written  as  the  product  of  the  individual 
translational,  rotational,  vibrational,  electronic  and  zero  point  contributions: 


Qs-mrs"Qf&sQ(s 


(22) 


The  translational  partition  function  reads: 


Ql  = 


2nMskT 

nJ2  , 


V 


(23) 


The  electronic  partition  function  is  computed  by  summing  up  the  contribution  from  each  excited  state: 


Q‘s  =  2^5 .*  eXP 

m 

The  zero  point  partition  function  is  given  by: 


f  e:  A 

S,m 


(24) 


Qs  =exp 


A°^ 

T 

V  1  J 


(25) 
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The  rotational  partition  function  for  rigid  rotation  of  a  linear  molecule  (N2,  H2,  02  and  OH)  is: 


Qs 


1  T 


The  rotational  partition  functions  for  nonlinear  molecules  (H20,  H02,  H202)  is  given  by: 


(26) 


QT 


—  I  ^ 

apr^r^r’3 


(27) 


The  vibrational  partition  function  of  the  harmonic  oscillator  is  computed  from  the  individual  contributions 
of  each  mode: 


qt= n 


i 


f 

l-exp 

V 


rivib,v  \ 

V S 


T 


J 


(28) 


The  derivation  of  the  relationship  between  the  partition  function  Q  and  the  thermodynamic  properties  is 
straightforward  ([1],  [91])  and  not  repeated  here.  For  example,  the  internal  energy  is  given  by: 


resulting  in: 


es=RT 


5  In  Qs 
dT 


(29) 


Zs-tfexpH^/D 

IT)-\+R  ^gfexpC-^'/T)  (3 

e 

With  Krot  expressing  the  number  of  rotational  degrees  of  freedom  (=2  for  linear  molecules  and  3  for  the 
nonlinear  ones). 

The  specific  heat  at  constant  volume  of  species  S  is  given  by: 


e?  —  — RT  +  K  r, 


RT 


-  + 


*1 


ef- 


exp(#5 


vib,\ 


resulting  in: 


cv,s 


des 

dT 


(31) 


3  Krot 

cvs=-R  +  —^- 

v’s  2  2 


R  +  Rj] 


6T'V  /  2 T 


sinh(#5  v  /  2 T) 


~\~  Cj7 


cel  = 


R 


Xg;'(^;')2exp  (-oeel/T) 


R 


T2  ^gfexpi-df  /T)  T2 


X^^exp(-^/r) 

^gfexpC-^/T) 


\2 


V  e 


(32) 


Finally,  the  entropy  can  be  expressed  as: 
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f 


ss=R 


kQ 


1  +  ln— 

v  R  j 


+  RT 


<3  In  Q 
dT 


(33) 


The  energy,  entropy  or  specific  heat  of  the  mixture  is  obtained  by  the  mass  fraction  weighted  average: 


(34) 


4.2.4  Laminar  transport  coefficients 

The  laminar  viscosity  of  each  individual  species  is  calculated  using  the  curve  fits  given  by  Blottner  [26]: 


Ms  =  1 — 5-exp(Cs  )T{As  ht)+Bs)  (35) 

m 

The  laminar  viscosity  of  the  gas  mixture  is  subsequently  computed  using  the  mixture  rule  of  Wilke  [95] 
which  is  suitable  for  reacting  gas  mixtures  without  ionisation  [64],  [6]: 


;  where 


n/s 


\  Ms ' 
1  +  — - 
M 


S  J 


1+ 


Ms 


J 


i 

^  f 


M, 


\MS  j 


(36) 


The  thermal  conductivity  of  the  individual  species  is  computed  by  the  Eucken  correction  utilizing  the 
modified  version  of  Hirschfelder  [16],  [36]: 


K S  —  Ms 


V  v  ,  (cv)T+(cv)T+(cvys 

2^v)s  Sc 


(37) 


The  heat  conductivity  of  the  mixture  is  computed  using  the  mixture  rule  of  Heming  and  Zipperer  [35]: 


k  = 


y  nsKs 


(38) 


4.2.5  Turbulence  modelling  and  Favre  averaged  conservation  equations 

Turbulent  boundary  and  shear  layers  are  present  in  flows  which  are  characterized  by  large  Reynolds 
numbers.  The  modelling  of  turbulent  flows  is  a  crucial  part  of  the  CFD  analysis  of  many  reacting  high 
enthalpy  flow  problems  of  practical  interest.  Because  of  the  complex  modelling  needed  to  account  for  the 
chemistry  and  thermodynamics,  the  application  of  advanced  turbulence  modelling  strategies  such  as  large 
or  detached  eddy  simulation  result  in  very  time  consuming  computations  which  are  currently  not  feasible 
for  design  applications  or  comprehensive  parametric  studies.  One  and  two  equation  eddy  viscosity  models 
combined  with  the  Favre  averaged  conservation  equations  are  well  established  and  represent  a  method  to 
model  turbulent  flows  in  a  robust  and  computationally  efficient  way.  The  solution  of  the  resulting  set  of 
equations  represents  a  time  average  of  the  flow  quantities.  Turbulent  fluctuations  are  not  resolved. 

The  averaged  equations  as  used  in  the  TAU  code  can  be  derived  from  the  laminar  or  instantaneous  set  of 
equations  (10)  to  (16)  by  applying  Reynolds  decomposition  to  the  total  and  partial  densities  and  Favre 
decomposition  to  the  remaining  flow  variables: 
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1  b+At 

p  =  p  +  p',p  =  —  [  p(t)dt,  (39) 

At  J 

h 

where  p  is  the  Reynolds  average  in  the  time  domain  of  the  density,  p ,  and  p'  is  the  fluctuating 
component  of  p. 

The  Favre  decomposition  of  the  remaining  flow  quantities  $  -  ( u ,  E ,  h,  T)  uses  a  density  weighted  average 
in  the  time  domain: 


--  t0+At 

</>  =  </>  +  f ,  &  =  —  |  p(t)</>(t)dt  (40) 

PAt  t 

The  integration  time,  At,  must  be  significantly  larger  than  the  largest  characteristic  time  of  the  considered 
fluctuating  turbulent  structures  [94].  If  this  is  not  the  case,  the  solution  of  the  averaged  equations  is  not 
unique  and  becomes  dependant  on  At  [23].  The  treatment  of  unsteady  flows  (URANS)  therefore  requires 
that  the  characteristic  time  scale  of  the  bulk  flow  is  sufficiently  larger  than  the  characteristic  time  scales  of 
the  turbulent  fluctuations. 


The  substitution  of  the  decompositions  (39)  and  (40)  in  the  laminar  or  instantaneous  set  of  conservation 
equations  (10)  to  (16),  results  in  the  Favre  averaged  conservation  equations  for  compressible,  turbulent 

flows.  The  turbulent  set  of  equations  then  contains  first  (e.g.  ui ),  second  (e.g.  pii-u" )  and  third  (e.g. 

puyy; )  moments  of  the  decomposed  fluctuating  flow  quantities.  The  derivation  of  transport  equations 

for  the  additional  terms  and  their  introduction  into  the  conservation  equations  would  result  in  the 
generation  of  further  unknowns  [94].  Therefore,  additional  turbulence  modelling  is  needed  to  close  the 
averaged  system  of  conservation  equations.  The  basis  of  Boussinesq  turbulence  models  is  the  introduction 
of  an  additional  eddy  viscosity,  pt,  which  is  computed  from  auxiliary  model  equations.  This  eddy  viscosity 
is  used  to  model  the  a-priory  unknown  higher  moments  of  the  decomposed  flow  quantities. 


Species  conservation  equations:  Turbulence  increases  the  species  diffusion  velocities  which  are 
modelled  by  simply  adding  the  turbulent  contribution  to  the  laminar  diffusion  coefficient: 


Ps—s  ~  P^^L 


f  Ps^ 


P__ 

Sc 


Pt 

Sc 


r  Ps 


(41) 


t  y 


Momentum  equations:  The  second  moments  which  are  introduced  by  the  averaging  procedure 
( puy' )  act  as  an  additional  turbulent  stress  tensor  (Reynolds  stress  tensor)  which  is  modelled 

using  the  Boussinesq  hypothesis  ( k  is  the  turbulent  energy,  the  last  term  represents  an  increase  of 
the  effective  static  pressure  due  to  the  presence  of  turbulent  fluctuations  and  is  neglected  for  the 
application  of  the  Spalart-Allmaras  one  equation  turbulence  model): 


-f  m  (42) 

•  Energy  equation:  The  definition  of  the  total  energy  changes  due  to  the  additional  contribution 
resulting  from  the  turbulence  to:  E  =  e  +  0.5u  •  u  +  k  .  The  turbulent  species  diffusion  is  modeled 
according  to  equation  (41)  and  the  turbulent  heat  conductivity  is  computed  assuming  a  ratio  of 
laminar  to  turbulent  Prandtl  numbers  of  Pr/PrT  =  0.8  using: 


PT  —  pT 


W  +  (  Vh7  )7 
yL 


-|(Yr«)/ 
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kt 


Pr  uT 

- !-j-k 

Pr7  jU 


(43) 


Introducing  the  additional  modeled  turbulent  contributions  from  equations  (41)  to  (43)  into  the  laminar  set 
of  equations  (10)  to  (16),  leads  to  the  modeled  Favre  averaged  set  of  conservation  equations. 


The  structure  of  the  Navier  Stokes  equations  in  their  integral  form  remains  unchanged: 


d_ 

dt 


\lMV  +  \l=undS=\l=sndS  +  \QdV 

VS  S  V 


The  vector  of  conservative  variables  now  contains  the  averaged  quantities: 


(44) 


y_=[ps,pu  ,pe) 

The  matrix  of  the  inviscid  fluxes  remains  unchanged: 


(45) 


Ps<± 

"  0 

II 

M 

puu 

+ 

Pi 

pEu 

- 1 

E-h 

*31 
I^H 
_ 1 

(46) 


The  matrix  of  the  viscous  fluxes  now  contains  the  additional  Reynolds  stresses  and  contributions  from  the 
turbulent  diffusion  and  heat  conduction.  Note  that  for  practical  computations  the  last  term  of  the 
approximation  of  the  Reynolds  stress  tensor  (42)  is  added  to  the  static  pressure:  peff  =p  +  2/3pk  and  is 

therefore  treated  as  an  inviscid  flux  contribution. 


y-iNS 

F  = 


f  \ 

SC  SCr 


t  y 

P  +  Pr 


p 


(^  K  +  Kj  ^  V  T  + 


P,  +  _Pt 


Sc  Sc 


t  y 


(47) 


Here  crk  is  a  constant  of  the  applied  two-equation  turbulence  model  and  k  denotes  the  turbulent  kinetic 
energy.  The  source  term  vector  now  contains  the  averaged  chemical  source  term  cos .  This  can  be 
estimated  e.g.  using  an  assumed  PDF  approach. 


Q  = 


(0S 

0 

0 


(48) 


The  eddy  viscosity,  //T,  can  be  computed  using  a  wide  variety  of  one-  and  two  equation  turbulence  models. 


Many  eddy  viscosity  turbulence  models  are  known  to  over  predict  the  spreading  rate  of  compressible  shear 
layers  above  convective  Mach  numbers  of  about  0.5  [94].  The  convective  Mach  number,  Mc,  is  related  to 
the  speed  at  which  large  turbulent  structures  are  transported  in  a  shear  layer  between  two  parallel  flows 
with  velocities  ux  and  u2  and  the  speeds  of  sound  rq  and  a2: 
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Mc  _  Ml - «2_  (49) 

Cl\  +  0,2 

Several  shear  layer  compressibility  corrections  exist  for  two-equation  eddy  viscosity  models  to  remedy 
this  problem.  Sarkar  et  al.  [76]  and  Zeeman  [97]  postulated  corrections  which  are  based  on  the  effects  of 
additional  dilatation  dissipation  and  eddy  shocklets  in  highly  compressible  flow.  Later,  Vreman  et.  al  [92] 
pointed  out  that  these  effects  are  not  significant  even  in  highly  compressible  mixing  layers.  Wilcox  [94] 
constructed  a  model  for  shear  layer  compressibility  correction  based  on  an  assessment  of  existing  work 
and  empirical  considerations  which  is  designed  to  minimize  the  adverse  effects  of  the  shear  layer 
correction  on  the  prediction  of  boundary  layer  flows. 


4.2.6  Modelling  of  flows  in  thermal  non-equilibrium 

The  numerical  analysis  of  e.g.  strongly  expanding  flows  in  nozzles  of  hypersonic  high  enthalpy  test 
facilities  requires  the  modeling  of  thermal  non-equilibrium  and  the  resulting  freezing  of  the  vibrational 
relaxation  process  of  the  molecules  in  the  flow  field.  This  is  achieved  in  TAU  by  considering  an  additional 
transport  equation  for  the  vibrational  energy  of  each  molecule  to  be  considered  in  thermal  non-equilibrium 
(N2,  02,  NO  for  air  flows).  In  accordance  with  the  nomenclature  of  equation  (44),  the  additional  fluxes  and 
source  terms  can  be  written  as: 


FEu  =[pefuT],  F A 


r 

\ 

Kv‘bVTTvib  + 

M  Mr 

£vibyT  Ps_ 

—  J 

KSc 

SCrj,  j 

b  - 

P  \ 

G=[nf 


(50) 


Further,  the  viscous  transport  of  vibrational  energy  by  conduction  has  to  be  considered  in  the  computation 
of  the  viscous  fluxes  in  the  transport  equation  of  the  total  energy  (47): 


^  total  energy 


(  K  +  Kj  )  V  T  + 


M  ,  Mr 
^  Sc  ScT  j 


Z4vr4+fe+^kf  +E*f  vTrsib 


(51) 


The  mixture  specific  vibrational  energy  of  a  species,  esvlb,  can  be  computed  from  the  specific  vibrational 
energy  using: 


pef=pssf>  (52) 

The  vibrational  temperature,  Tsvlb,  is  related  to  the  vibrational  energy  e.g.  through  the  set  of  equations  (22) 
to  (33)  defining  the  species  partition  functions.  The  vibrational  heat  conductivity,  Ksvlb,  is  computed  using 
the  modified  Eucken  correction  and  applying  turbulence  and  mixing  rule  factors  similar  those  given  in 
equations  (53)  and  (54): 


U-Vib  - 
KS  ~ 


l  +  ZLit 

Prr  p 


ns  (cv  ), 


vib 


0s  Sc 


-Ms 


(55) 


The  source  term  of  the  vibrational  energy  equations  consists  of  the  chemical  contribution  (ss  0)s)  and  the 
contribution  from  the  vibration-translation  relaxation  (I2svlb,t). 


The  vibration-translation  relaxation  term ,nsvlb,t,  is  computed  using  the  Landau-Teller  approach  with  the 
equilibrium  energy  set  to  the  vibrational  energy  at  the  translational  temperature,  T: 
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(56) 


The  translation- vibration  relaxation  times  of  different  pairs  of  species  s,S  are  computed  using  the  empirical 
relation  of  Millikan  and  White: 
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(57) 


4.2.7  Modelling  of  weakly  ionized  flows 

The  assumptions  of  an  independent  harmonic  oscillator  and  rigid  rotator  which  are  made  for  the 
calculation  of  the  partition  function  of  molecules  are  often  questionable  at  large  temperature  levels. 
Improved  modeling  of  the  thermodynamic  properties  of  molecular  species  can  be  achieved  by  considering 
a  full  coupling  of  the  rotational,  vibrational  and  electronic  excitation  as  well  as  anharmonicity  effects.  A 
detailed  description  of  this  approach  is  given  e.g.  in  reference  [6].  Exemplary,  the  internal  partition 
function  of  a  diatomic  molecule  can  be  expressed  as: 


Qint  =TjSe  (-0JT) 

e 


jr _ l 

(j0r  l-exp  (-0JT) 


x 


1  + 


eB  +  ol 


3  T  15  T1 


-  + 


2PJ  ae _ 

~ba +s„(i-exp(0v/:r))  + 


XqBJT 

(exp(yT)-l) 


(58) 


The  first  terms  in  equation  (58)  are  the  harmonic  contributions  of  the  electronic,  rotational  and  vibrational 
partition  functions  with  the  characteristic  temperatures  9^  and  9W.  The  last  term  models  the  coupling 
between  the  internal  energy  modes  and  involves  the  additional  coupling  constants  ae,  x0  and  B0.  If 
polyatomic  molecules  are  considered,  the  rotational  term  has  to  be  modified  and  further  vibrational  terms 
for  each  additional  vibrational  mode  have  to  be  added  (see  e.g.  Gurvich  [27]). 


Concerning  the  computation  of  transport  coefficients,  the  Yos  mixture  rule  is  generally  considered  to  be 
more  accurate  than  the  Wilke  rule  for  high  temperatures  and  the  presence  of  ionizes  species.  It  is  an 
approximation  of  the  Chapman-Enskog  expansion  for  the  transport  coefficient  which  neglects  off-diagonal 
terms  (details  and  derivation  in  ref.  [6]).  Different  relations  exist  for  the  viscosity  and  the  translational  and 
internal  thermal  conductivities.  The  computation  procedure  for  the  viscosity  is  repeated  here  from  [6]  as 
an  illustrative  example  (. x  are  the  species  mole  fractions): 
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(60) 

(61) 


The  binary  collision  integrals  in  the  set  of  equations  (59)  to  (61)  are  obtained  from  curve  fits.  For  neutral 
particles  the  curve  fit  function  is  given  in  equation  (62): 


InQ*  =D  +  (Aln2r  +  filnr  +  C)lnr  (62) 

The  coefficients  A,  B ,  C  and  D  are  needed  for  each  pair  of  species  (ij).  Hence,  the  procurement  of  the 
required  input  data  (especially  for  a  large  number  of  species)  is  significantly  more  complex  than  for  the 
Wilke  /  Blottner  approach. 

The  application  of  Fick’s  law  to  high  temperature  gas  mixtures  involving  free  electrons  and  ionized 
species  leads  often  to  inaccurate  results  for  the  diffusive  mass  fluxes.  Ramshaw  and  Chang  [69],  [70] 
developed  an  appropriate  ambipolar  diffusion  model  for  ionized  mixtures.  It  was  demonstrated  in  [25]  that 
for  re-entry  applications  this  “Self  Consistent  Effective  Binary  Diffusion  Method”  (SCEBDM)  leads  to 
similar  results  as  the  solution  of  the  full  Stefan-Maxwell  equations. 

The  diffusive  mass  flux  of  species  S  is  computed  according  to: 
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where  c  is  the  total  molar  concentration  of  the  mixture  and  q  are  the  species  charges.  The  effective 
diffusion  coefficient  Djm  for  species  j  is  computed  from  the  set  of  binary  diffusion  coefficients  Djt  using: 


i  species 


D,.=(!  -A/Xf- 

/  j=\,j*i  Uij 


(64) 


4.2.8  Numerical  solution  procedure,  applied  spatial  and  temporal  discretization  schemes 

The  spatial  discretization  used  in  the  TAU  code  is  based  on  a  cell-vertex  scheme.  The  primary 
computational  grid  consists  of  tetrahedra,  prisms,  hexahedra  and  pyramids.  The  flow  quantities  are  stored 
on  the  vertices  of  the  primary  grid.  For  the  finite  volume  discretization,  a  dual  grid  cell  is  constructed 
around  each  vertex  of  the  primary  grid.  A  schematic  of  this  arrangement  is  shown  in  Figure  24. 
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Figure  24:  Schematic  of  the  dual  cell  structure  applied  for  the  spatial  discretization  in  TAU 

The  dual  cell  with  volume,  Vh  and  surface,  Su  which  is  constructed  around  point  PL  in  Figure  24  is  marked 
in  grey.  The  faces  between  two  adjacent  volumes  are  characterized  by  their  normal  vectors.  They  are 
grouped  such  that  one  face  is  always  associated  to  two  adjacent  points  in  the  primary  grid.  For  the 
example  in  Figure  24,  the  normal  vector,  n ,  of  the  face  which  connects  the  dual  cells  around  the  points  PL 
and  PR  is  computed  from  its  components  using: 


n  =  nx  +  n2  (65) 

The  boundary  of  the  dual  grid  cells  is  constructed  by  connecting  the  centres  of  the  primary  grid  cells  (e.g. 
point  B)  with  the  centres  of  the  primary  grid  faces  (e.g.  point  C).  The  centres  of  the  primary  grid  cells  are 
derived  from  a  weighted  average  of  the  vertices.  For  the  two-dimensional  example  in  Figure  24  the 
coordinates  of  point  B  are  computed  using: 


B  = 


vv,  PL  +w2A  +  w3Pr 


(66) 


W]  +  w2  +  w3 

The  weighting  functions  represent  the  sum  of  the  squared  face  areas  or  side  lengths  which  are  adjacent  to  a 
vertex  point,  i.e.: 


^=(pLp«)2 +(pla)2  (67) 

Subsequently,  the  discretized  form  of  the  set  of  conservation  equations  (44)  can  be  written  as: 

—  \udV:=V—U,\FEu/NSndS:=YFEu/NSni,  f QdV:  =  QV  (68) 

r)t  *  dt  ^  J  —  — 

v  dV  1  V 

The  discretized  conservation  equations  for  one  control  volume  with  the  face  segments,  /,  in  the  final  form 
introducing  the  residual,  R ,  are: 
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dr 


-U  =  -R 


R  =  — 
-  V 


'Ll** 

V  i  i 


(69) 


The  actual  computation  of  the  inviscid  flux  vector,  n{ ,  across  a  cell  face,  i,  is  done  using  a  variety  of 
upwind  or  central  schemes. 


To  increase  the  accuracy  of  the  upwind  discretization  to  second  order  in  space,  the  flow  quantities  at  the 
left  and  right  side  of  each  face,  (j)L  and  (j)R,  are  extrapolated  assuming  a  piecewise  linear  distribution  using: 


^=^+Uw<t>L{xR-xL) 

<Pfr  -fRode +\mR(xL-xR) 


(70) 


The  vectors  xL  and  xR  represent  the  coordinates  of  the  dual  cell  centres  of  two  adjacent  volumes  (points 
PL  and  PR  in  Figure  24).  The  limiter  function,  L,  is  introduced  to  prevent  overshoots  of  the  reconstructed 
quantities  due  to  the  presence  of  discontinuities.  In  the  unstructured  regions  of  the  grid  the  Barth- Jespersen 
limiter  [3]  is  used. 


The  gradients  of  the  flow  quantities  which  are  needed  for  the  piecewise  linear  reconstruction  to  achieve 
second  order  accuracy  in  space  for  the  upwind  discretization  and  for  the  computation  of  the  viscous  fluxes 
are  computed  using  a  least-square  fitting  approach  as  described  by  Anderson  [2]  and  Haselbacher  [34]. 
The  flow  quantities  at  the  grid  points  surrounding  the  cell  for  which  the  gradients  are  to  be  computed  are 
approximated  using  a  second  order  Taylor  expansion  and  the  gradients  in  this  expansion  are  chosen  such 
that  the  resulting  estimate  provides  an  optimum  least-square  fit  to  the  actual  computed  cell  averages.  In 
structured  regions  with  aligned  cells  of  large  aspect  ratio  (e.g.  meshes  close  to  viscous  walls  consisting  of 
hexahedra  and  prisms)  the  cell  boundaries  form  grid  lines  with  minor  change  of  direction.  Those  lines  are 
identified  and  the  gradients  along  these  grid  lines  are  computed  using  a  second  order  finite  difference 
approach.  The  van-Leer  limiter  [88]  is  used  for  the  correction  of  these  gradients. 


The  viscous  flux  vector,  F_^s  ,  across  a  cell  face,  /,  is  computed  according  to  the  equations  given  above 

(e.g.  equation  (47)  for  the  viscous  fluxes  in  the  continuity,  momentum  and  total  energy  equation).  The 
flow  quantities,  ((),  and  the  gradients  of  the  flow  quantities,  V(|),  at  the  cell  faces  which  are  needed  in  these 
equations  (see  also  Figure  24)  are  computed  from  the  averaged  quantities  of  the  dual  grid  cells  adjacent  to 
the  considered  face  using  a  central  scheme: 


<t>face=\(<t>,+(t>R)  ,V^=|(V&  +  V&)  (71) 

The  accuracy  of  the  gradient  reconstruction  for  the  computation  of  the  viscous  fluxes  is  improved  by 
applying  the  following  correction  procedure: 


V</>- 


face, corrected 


=  v</>face  + 


^LLR 

V 

1 

¥ 

5 

!  J 

\^*LR  | 

(72) 


with  Axlr  denoting  the  spatial  difference  vector  between  the  left  and  right  dual  cell  centres  belonging  to 
the  face  at  which  the  gradient  needs  to  be  evaluated.  This  procedure  corrects  the  component  of  the  vector 
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tangential  to  the  spatial  difference  vector  to  the  central  difference  value  while  preserving  the  normal 
components. 

To  enhance  the  convergence  rate  of  the  iterative  solution  of  equation  (69),  explicit  smoothing  of  the 
residuals  or  multigrid  schemes  can  be  applied. 

Temporal  discretization  of  equation  (69)  is  achieved  by  applying  multi-stage  Runge  Kutta  schemes  as 
described  e.g.  by  Jameson  et  al.  [40]  or  an  implicit  approximate  factorization  scheme  (Lower-Upper 
Symmetric  Gauss-Seidel).  For  steady  state  problems,  equation  (69)  reduces  to  R= 0  because 
dU_/dt  vanishes.  The  iterative  solution  is  performed  in  a  pseudo  time,  t*\ 


d  U 
d  t 


+  R  =  0 


(73) 


When  explicit  Runge-Kutta  time  integration  is  applied,  the  computation  of  the  chemical  source  term 
vector  is  done  in  a  point  implicit  manner  as  described  by  Sheffer  [80].  The  vector  of  the  chemical  sources 
is  linearized  to  obtain  an  estimate  for  the  new  time  step  (n+ 1): 


Q 


n+ 1 


:G»+^L  — l/Af* 

-  8U  dt~ 


Substitution  into  the  governing  equation  yields: 


(74) 


—  f  UdV  = 

8t{- 

which  is  evaluated  entirely  at  the  current  time  level  and  is  thus  fully  explicit.  This  treatment  necessitates 
the  inversion  of  the  source  term  Jacobian  matrix  with  dimension  NxN  where  N  is  the  number  of  species. 
The  dependency  of  the  source  term  on  the  temperature  is  neglected  and  no  stability  penalties  from  this 
approximation  are  observed. 

Time  accurate  computations  are  performed  using  a  dual-time  stepping  technique  as  described  by  Jameson 
[39].  The  computational  time  span  [0;  t |  is  subdivided  into  intervals  At  (“physical  time  step”).  The 
application  of  second  order  backward  differencing  in  time  to  the  vector  of  conservative  variables  in 
equation  (69)  yields: 


.  8Qn 

I  At  ~ 

-t 

f  FEundS  +  [FNSndS  +  f  OdV 

8U 

1 

II 

i 

II 

1 

II 

1 _ 

(75) 


3  U  2  U  U 

— - —  +  = —  =  -R 

2A  t  At  2A  t 

If  the  solutions  at  the  time  steps  n  and  n- 1  are  known,  equation  (76)  can  be  rearranged  to: 


(76) 


d  Un+l 


i  * 

dt 


+  Ra 


=  0 


j^dual 


=  Rn+l 


3(7  "+1 

H - 

2A  t 


2(7" 

At 


H - 

2A  t 


(77) 


The  standard  solution  procedure  for  steady  state  problems  including  all  means  for  convergence 
acceleration  can  then  be  used  to  solve  equation  (77)  in  pseudo  time,  t*,  to  advance  the  solution  to  the 
physical  time  step  n+ 1. 
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The  computational  time  step,  A t*,  is  computed  locally  for  each  cell  (“local  time  stepping”)  using  the 
minimum  of  the  convective  and  viscous  time  steps: 


A t  =min(A£c ,  A^v)  (78) 

Both  time  steps  are  determined  by  adding  the  contributions  from  the  respective  convective  and  viscous 
Eigenvalues  from  all  cell  faces,  i. 


&1‘  =CFL  2V=!—  ,  A f  =CFL-2V2  —  (79) 

2>;  4  £4 

i  i 

The  contribution  from  the  maximum  convective  Eigenvalue  evaluated  at  each  cell  face,  /,  is: 


K  =“/«/+«/  ki  |  (80) 

The  contribution  from  the  viscous  Eigenvalues  is  estimated  as: 


4=(4+4)hf  (8i) 

The  Eigenvalues  and  Avl  correspond  to  the  contributions  from  viscous  momentum  transport  and  the 
heat  conductivity: 


r  _  4 

4-- 


JU  +  /uT 

V  P 


Q  i  _  P 

’  A/2  —  — 

CvP 


f  \ 

M  +  Mr 


Pr  Pr7 


t  7 


(82) 


The  contribution  from  species  diffusion  to  the  viscous  time  step  is  neglected  because  it  is  usually  (for 
moderate  Schmidt  numbers  greater  than  0.6)  smaller  than  the  contribution  from  heat  conductivity: 


J!_+Jh_<xv2 

pSc  pScT 


(83) 


5.0  EXAMPLES  OF  COMBINED  EXPERIMENTAL  AND  NUMERICAL  FLOW 
FIELD  INVESTIGATIONS 

5.1  Calibration  of  HEG  Free  Stream  Flow 

The  high  enthalpy,  high  pressure  flows  generated  in  free  piston  driven  shock  tunnels  are  well  suited  for  the 
investigation  of  chemically  reacting  flows  occurring  during  the  (re-)entry  flight  of  space  vehicles  or  in 
combustors  of  hypersonic  airbreathing  propulsion  system.  The  proper  utilisation  of  such  ground  based 
impulse  facilities  requires  a  detailed  knowledge  of  the  test  section  flow.  However,  due  to  the  flow 
conditions  generated  in  high  enthalpy  tunnels,  the  calibration  process  is  much  more  complex  than  in  cold 
hypersonic  facilities.  One  reason  for  this  is  that  the  behaviour  of  the  test  gas  deviates  from  that  of  a 
calorically  perfect  gas.  With  increasing  total  enthalpy,  the  excitation  of  internal  degrees  of  freedom  such 
as  the  vibrational  excitation  of  molecules  and  subsequently  chemical  reactions  must  be  considered. 
Additional  difficulties  arise  from  the  fact  that  some  parts  of  the  HEG  facility,  such  as  the  nozzle  reservoir, 
are  not  accessible  with  modern  optical  measurement  techniques  due  to  limitations  caused  by  their 
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structural  design.  These  result  from  temperatures  of  up  to  10000  K  and  pressures  of  up  to  2000  bar  which 
are  reached  at  the  highest  total  enthalpy  and  total  pressure  conditions.  Therefore,  the  evaluation  of  the  free 
stream  conditions  in  the  test  section  of  the  HEG  free  piston  driven  shock  tunnel  is  performed  by  numerical 
analysis  using  a  suitable  set  of  measured  input  and  reference  parameters.  The  general  iterative  procedure 
established  to  determine  HEG  free  stream  conditions  is  outlined  in  Figure  25.  The  numerical 
determination  of  the  HEG  free  stream  consists  of  two  steps  which  require  a  suitable  set  of  input 
parameters.  First,  the  nozzle  reservoir  temperature  is  computed  with  a  one-dimensional  analysis  of  the 
shock  tube.  The  relevant  input  parameters  are  the  measured  values  of  the  initial  shock  tube  filling  pressure 
and  temperature,  the  shock  speed  and  the  nozzle  reservoir  pressure.  Based  on  these  nozzle  reservoir 
conditions,  the  free  stream  is  subsequently  determined  by  numerical  simulation  of  the  nozzle  flow. 


Shock  tube  measurements 

•  Reservoir  pressure 

•  Initial  pressure,  temperature 
in  shock  tube 

•  Shock  speed 


Nozzle  measurents 

•  Wall  heat  flux  (Boundary 
layer  transition  point) 

•  Wall  pressure 


Free  stream  measurements 

•  Static  pressure 

•  Pitot  pressure 

•  Stagnation  point  heat 
transfer  on  a  spherical  probe 


Shock  tube  computation 


I 


Reservoir  conditions 


1 


Nozzle  flow  computation 


I 


Free  stream  conditions 


Cross  -  check 


& 

T3 

O 


E 


Figure  25:  Combined  experimental  and  numerical  iterative  procedure  to  determine  HEG  free 

stream  conditions 

The  close  cooperation  between  experiment  and  CFD  is  an  interactive  process  which  helps  to  validate  new 
measurement  techniques  as  well  as  the  physico-chemical  model  used  in  the  CFD  tools  such  providing  step 
by  step  an  improved  understanding  of  the  facility  performance.  In  free  piston  shock  tunnels  the  available 
test  time  is  in  the  order  of  milliseconds.  In  order  to  be  able  to  select  the  useful  test  time,  the  starting 
process  of  the  nozzle  flow  and  the  flow  past  the  model  and  the  arrival  of  either  the  expansion  wave  which 
reflects  off  the  piston  face  or  the  helium  driver  gas  in  the  test  section  must  be  known.  In  HEG  the 
determination  of  the  test  time  window  -  in  particular  for  the  high  enthalpy  conditions  -  is  based  on  the 
evaluation  of  the  measured  time  resolved  nozzle  reservoir  pressure  as  well  as  the  Pitot  pressure  and  the 
static  pressure  in  the  test  section. 

5.1.1  High  Enthalpy  Operating  Conditions  I-IV 

For  the  high  enthalpy  (12-23MJ/kg)  operating  conditions  of  HEG,  nozzle  2  is  utilised.  This  conical 
nozzle  has  an  expansion  half-angle  of  approximately  6.5°,  a  total  length  of  3.75  m,  a  throat  radius  of 
0.01 1  m  and  an  exit  radius  of  0.44  m.  The  area  ratio  of  the  nozzle  is  1600.  The  throat  section  was  designed 
using  fourth  and  third  order  polynomials  in  order  to  obtain  a  smooth  nozzle  contour,  i.e.  a  contour 
whereby  the  radius  as  a  function  of  the  longitudinal  coordinate  has  continuous  first  and  second  derivatives. 

The  evaluation  of  the  heat  transfer  measurements  along  the  nozzle  wall  indicates  that  a  transition  region 
exists  between  0.5  -  1.0  m  downstream  of  the  nozzle  throat  [30].  In  order  to  investigate  the  influence  of 
the  location  of  the  transition  point  on  the  shape  of  the  computed  Pitot  pressure  profiles  at  the  nozzle  exit, 
the  normalised  measured  and  computed  Pitot  pressure  profiles  are  plotted  for  operating  condition  I  in 
Figure  26.  It  should  be  noted  that  these  numerical  results  were  obtained  using  the  CEVCATS-N  code.  The 
measured  Pitot  pressure  was  evaluated  in  two  different  test  time  windows  and  no  significant  change  of  the 
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shape  was  found.  The  computed  Pitot  pressure  profile  shapes  depend  strongly  on  the  modelling  of  the 
boundary  layer  flow.  A  fully  laminar  and  two  turbulent  computations  with  different  locations  of  the 
laminar-turbulent  transition  point  were  performed.  For  condition  I,  the  best  agreement  between  the 
measured  and  computed  shape  could  be  obtained  when  assuming  that  transition  occurs  at  the  downstream 
end  of  the  transition  region  found  experimentally. 


r[m] 


Figure  26:  Comparison  between  measured  and  computed  Pitot  pressure  profiles  showing  the 
influence  of  the  location  of  the  laminar  to  turbulent  transition  point;  HEG  condition  I 


The  numerical  sensitivity  analyses  included  computations  assuming  a  chemical  /  thermal  nonequilibrium, 
chemical  nonequilibrium  /  thermal  equilibrium  and  chemical  /  thermal  equilibrium  nozzel  expansion.  The 
chemical  /  thermal  nonequilibrium  computations  were  additionally  performed  applying  different  sets  of 
chemical  reaction  rates.  The  chemical  /  thermal  nonequilibrium  computation  using  the  reaction  rate  data 
given  by  Park  [65]  was  chosen  as  the  reference  computation.  Compared  to  this  computation,  the  results 
obtained  using  the  reaction  rates  published  by  Dunn  and  Kang  [18]  and  Gupta  et  al.  [26]  tend  to  slow 
down  the  chemical  relaxation  process  during  the  expansion.  This  results  in  a  reduction  of  the  exit 
translational  temperature  by  5%  and  23%,  respectively  (Figure  27). 


Figure  27:  Sensitivity  analysis  concerning  the  computed  temperature  distribution  along  the 
nozzle  centerline,  HEG  condition  I;  transition  location  at  x=0.75m 

While  the  Pitot  pressure  remains  approximately  constant,  the  static  pressure  decreases  by  4%  and  16%  and 
the  density  increases  by  0.5%  and  3.5%.  When  comparing  the  chemical  /  thermal  nonequilibrium  and  the 
chemical  nonequilibrium  /  thermal  equilibrium  computations,  the  exit  temperature  is  increased  by  22%, 
the  Pitot  pressure  is  reduced  by  4%,  the  static  pressure  is  increased  by  17%  and  the  density  is  reduced  by 
5%.  This  trend  is  further  intensified  when  comparing  the  chemical  /  thermal  nonequilibrium  computation 
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with  the  one  assuming  an  expansion  in  chemical  /  thermal  equilibrium.  Here  the  temperature  is  increased 
by  160%,  the  Pitot  pressure  is  reduced  by  17%,  the  static  pressure  is  increased  by  83%  and  the  density  is 
reduced  by  20%. 


Figure  28:  Temporal  development  of  the  reservoir  pressure,  static  pressure  and  Pitot  pressure  in 
the  test  section  for  HEG  condition  I  (left)  and  condition  III  (right) 


Figure  29:  Comparison  of  measured  and  computed  Pitot  pressure  profiles  at  the  nozzle  exit.  The 
nozzle  flow  of  pure  air  as  well  as  of  air  plus  5%  and  10%  mass  fraction  helium  was  computed; 

HEG  condition  I 


The  temporal  development  of  the  reservoir  pressure,  the  Pitot  pressure  and  the  static  pressure  in  the  test 
section  for  operating  conditions  I  and  III  are  plotted  in  Figure  28.  In  order  to  allow  a  direct  comparison  of 
the  temporal  development  of  these  three  quantities,  the  reservoir  pressure  and  the  static  pressure  were 
scaled.  Therefore,  the  unit  of  the  ordinate  in  Figure  28  is  only  related  to  the  Pitot  pressure  signal.  The  time 
scale  on  the  abscissa  starts  when  the  incident  shock  in  the  shock  tube  reflects  off  the  end  wall.  The 
reservoir  pressure  signal  was  shifted  by  approximately  0.5  ms  for  conditions  I  and  0.7  ms  for 
conditions  III  in  order  to  account  for  the  delay  of  the  flow  start-up  in  the  test  section.  A  period  of  time  of 
approximately  1.3  -  1.5  ms  is  necessary  to  establish  a  steady  nozzle  flow.  After  the  nozzle  starting  process, 
the  Pitot  pressure  and  static  pressure  development  in  the  test  section  is  clearly  driven  by  the  reservoir 
pressure.  For  condition  I,  however,  the  static  pressure  decreases  relative  to  the  development  of  the 
reservoir  pressure  and  the  Pitot  pressure  after  about  1.8  ms.  This  drop  of  the  static  pressure  is  much  less 
pronounced  for  condition  III.  Here  it  starts  after  approximately  2.7  ms.  For  condition  IV  and  condition  II 
(not  shown  here)  this  effect  can  not  be  observed.  One  possible  explanation  of  this  behaviour  could  be  the 
arrival  of  helium  driver  gas  in  the  test  section. 
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In  order  to  study  the  effect  of  driver  gas  contamination  in  more  detail,  steady  state  nozzle  flow 
computations  were  performed  for  condition  I  using  a  test  gas  mixture  of  air  and  helium.  For  these 
computations  a  homogeneous  mixture  of  air  and  helium  is  used.  The  nozzle  reservoir  conditions  are 
determined  by  assuming  that  the  helium  has  the  same  enthalpy  as  at  the  end  of  the  compression  tube  at 
main  diaphragm  rupture.  This  results  in  approximately  the  same  reservoir  enthalpy  as  given  in  Table  1.  In 
Figure  29,  the  computed  Pitot  pressure  profiles  at  the  nozzle  exit  for  pure  air,  air  plus  5%  mass  fraction 
helium  and  air  plus  10%  mass  fraction  helium  are  compared  with  the  measured  profiles.  The  changes  of 
the  Pitot  pressure  profiles  due  to  the  addition  of  helium  are  small  when  compared  with  the  fluctuations  of 
the  measured  data.  In  contrast  to  the  Pitot  pressure,  the  static  pressure  is  significantly  reduced  when 
helium  is  added  to  the  test  gas  air.  The  averaged  free  stream  static  pressure  at  the  nozzle  exit  resulting 
from  nozzle  flow  computations  using  pure  air  and  air  plus  helium  is  plotted  in  Figure  30  together  with  the 
measured  temporal  development  of  the  static  pressure  for  different  condition  I  runs.  The  computed  static 
pressure  for  pure  air  agrees  well  with  the  measured  value  at  about  1.5  ms,  i.e.  after  the  nozzle  flow  starting 
process  has  finished.  After  a  period  of  time  of  approximately  0.2  -  0.3  ms  the  static  pressure  drops  and  the 
addition  of  helium  in  the  computations  has  the  same  effect.  A  concentration  of  5%  mass  fraction  of  helium 
reduces  the  static  pressure  by  approximately  20%. 

The  trends  provided  by  the  sensitivity  analyses,  described  above,  allow  to  choose  the  physico-chemical 
modelling  which  is  best  suited  to  reproduce  the  whole  set  of  measured  data.  In  general,  the  best  agreement 
with  the  experimental  data  was  obtained  assuming  the  nozzle  flow  to  be  in  chemical  nonequilibrium  and 
thermal  equilibrium.  This  assumption  is  further  supported  by  LIF  measurements  in  HEG  which  indicate 
that  the  vibrational  temperature  of  NO  is  in  equilibrium  with  the  translational  /  rotational  temperature 
(Wollenhaupt  et  al.  [96]). 


t  |ms| 


Figure  30:  Comparison  between  the  measured  temporal  development  of  static  pressure  and 
computed  averaged  static  pressure  in  the  test  section  core  flow  resulting  from  steady  state 
nozzle  flow  computations  using  pure  air  and  air  plus  5%  and  10%  mass  fraction  helium;  HEG 

condition  I 

In  addition  to  the  computation  of  test  section  free  stream  conditions  and  comparison  with  corresponding 
measurements,  the  investigation  of  the  flow  past  basic  aerodynamic  configurations  is  well  suited  to  study 
fundamental  aspects  of  high  enthalpy  flow  fields  as  well  as  to  validate  the  facility  performance, 
measurement  techniques  and  CFD  codes.  A  transverse  cylinder  with  a  diameter  of  Dcyi  =  90  mm  was 
chosen  in  order  to  further  investigate  the  influence  of  early  driver  gas  contamination  on  the  considered 
flow. 

For  HEG  condition  I,  computational  flow  field  analyses  were  first  performed.  In  order  to  carry  out  an 
assessment  of  the  physico-chemical  model  used  in  the  CFD  code,  computations  were  again  performed 
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assuming  a  chemical  nonequilibrium  /  thermal  equilibrium,  chemical  /  thermal  nonequilibrium  and  a 
frozen  flow.  For  the  chemical  nonequilibrium  computations,  the  reaction  rates  of  Park,  Dunn  &  Kang  and 
Gupta  were  again  used.  The  chemical  /  thermal  nonequilibrium  computation  showed  that  for  HEG 
condition  I  the  flow  reaches  a  thermal  equilibrium  state  a  short  distance  behind  the  bow  shock.  Compared 
to  a  chemical  nonequilibrium  /  thermal  equilibrium  computation,  the  resulting  difference  in  shock  stand¬ 
off  distance  amounts  to  less  than  2%  and  is,  therefore,  significantly  smaller  than  the  difference  obtained  by 
using  different  chemical  models.  The  computed  density  distributions  along  the  stagnation  streamline  in 
Figure  31  show  that  the  addition  of  10%  helium  also  causes  a  larger  shift  of  the  shock  stand  off  distance 
than  obtained  by  using  different  chemical  models.  Clearly,  when  helium  is  added,  the  flow  field 
characteristics  change  from  a  nonequilibrium  flow  with  chemical  relaxation  behind  the  bow  shock, 
causing  an  increase  of  density,  towards  a  frozen  flow.  As  a  reference,  the  density  distribution 
corresponding  to  a  frozen  shock  layer  flow  is  also  shown  in  Figure  3 1 . 

The  temporal  development  of  the  shock  stand-off  distance,  which  was  obtained  by  the  evaluation  of  time 
resolved  shadowgraphs,  is  plotted  in  Figure  32  together  with  the  development  of  the  static  pressure  and  the 
Pitot  pressure.  The  increase  of  the  shock  stand-off  distance  between  approximately  1.8  ms  and  2.2  ms 
coincides  with  the  drop  of  the  static  pressure.  It  is  interesting  to  note  that  for  t  >  2.8  ms  the  drop  of  the 
static  pressure  is  significantly  reduced  and  the  shock  stand-off  distance  remains  approximately  constant. 
This  would  imply  that  after  an  initial  amount  of  driver  gas  has  reached  the  test  section,  its  concentration 
stays  approximately  constant  and  is  not  continuously  increased  at  the  same  rate.  This  observation  is 
consistent  with  the  numerical  results  of  Chue  and  Eitelberg  [12],  who  computed  the  transient  flow  at  the 
end  of  the  HEG  shock  tube.  They  found  that  the  interaction  of  the  reflected  shock  with  the  contact  surface 
immediately  results  in  a  stream  of  driver  gas  towards  the  shock  tube  end  wall.  However,  further  leakage  of 
driver  gas  is  retarded  due  to  vortical  structures  generated  at  the  contact  surface.  They  conclude  that  the 
initial  jetting  of  the  driver  gas  into  the  test  gas  is  not  necessarily  sustained  and  that  the  subsequent  events 
could  be  more  associated  with  the  mixing  of  the  initial  leakage  gas  than  with  additional  contaminant. 


X/Rcyl 


Figure  31 :  Normalised  computed  density  distribution  along  the  stagnation  streamline  of  a 
cylinder  flow  in  HEG  condition  I;  comparison  of  different  sets  of  chemical  reaction  rates 


The  shock  stand  off  distance  and  the  shock  shape  resulting  from  a  schlieren  image  and  CFD  are  compared 
in  Figure  33.  The  measurement  was  taken  at  3.75  ms  after  shock  reflection,  i.e.  after  the  helium 
contamination  took  place.  Good  agreement  between  computed  and  measured  shock  stand-off  distance  and 
shock  shape  could  be  found  for  the  computation  using  an  air  /  10%  helium  mixture. 
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Figure  32:  Temporal  development  of  the  shock  stand-off  distance  in  comparison  with  the 
development  of  the  Pitot  pressure  and  the  static  pressure;  condition  I,  air 


Figure  33:  Comparison  of  measured  and  computed  shock  shapes,  HEG  operating  condition  I 

In  order  to  delay  the  early  driver  gas  contamination,  a  method  which  was  proposed  by  Sudani  et  al.  [87] 
was  implemented  in  HEG.  In  [87]  it  was  demonstrated  that  a  sleeve,  mounted  at  the  end  of  the  shock  tube 
of  the  T5  Hypervelocity  Shock  Tunnel  of  GALCIT  is  very  efficient  in  delaying  the  early  driver  gas 
contamination.  This  device  was  designed  in  such  a  way  that  it  captures  the  driver  gas  jet  along  the  wall 
upstream  of  the  shock  tube  end  wall  before  it  mixes  with  the  test  gas.  For  HEG  a  similar  device  was 
designed  (see  Figure  23)  based  on  computations  of  the  shock  reflection  process  in  HEG  by  Chue  and 
Eitelberg  [12].  As  shown  in  Figure  34,  a  sleeve  length  of  150  mm  was  chosen.  The  gap  between  the  shock 
tube  wall  and  the  sleeve  was  set  to  2  mm.  Six  110  mm  long  and  4  mm  thick  struts  which  are  evenly 
distributed  along  the  circumference  were  used  to  connect  the  sleeve  with  the  shock  tube  wall.  Due  to  the 
design  of  the  downstream  closure  of  HEG,  it  is  not  possible  to  suck  the  driver  gas  into  a  driver  gas 
removal  chamber  as  was  done  successfully  in  T5.  The  experiments  with  the  sleeve  were,  therefore, 
performed  by  only  using  the  volume  between  the  sleeve  and  the  shock  tube  plus  the  additional  volume 
above  the  shock  tube  end  plate  (see  Figure  34)  to  capture  the  helium. 
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Figure  34:  Sleeve  mounted  at  the  end  of  the  HEG  shock  tube;  measures  in  mm 


Figure  35:  Temporal  development  of  the  reservoir  pressure,  the  static  pressure  and  the  Pitot 
pressure  in  the  test  section  for  condition  I;  run  513  without  sleeve  (upper)  and  run  514  and  515 

with  sleeve  (middle  and  lower) 
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In  Figure  35,  the  temporal  development  of  the  reservoir  pressure  and  the  Pitot  pressure  and  static  pressure 
in  the  test  section  are  plotted  in  the  same  way  as  in  Figure  28.  The  oscillations  of  the  static  pressure  signal 
at  approximately  2.8  ms  is  caused  by  electromagnetic  interference  caused  by  a  laser  pulse  used  for  LIF 
measurements.  The  results  of  three  subsequent  runs  are  plotted.  The  upper  plot  refers  to  a  run  without 
sleeve  and  the  data  plotted  in  the  middle  and  lower  diagram  were  obtained  with  the  sleeve  mounted  in  the 
shock  tube.  It  is  obvious,  that  the  early  drop  of  the  static  pressure  relative  to  the  reservoir  and  Pitot 
pressure  does  not  occur  during  the  runs  with  sleeve.  Based  on  the  investigations  discussed  in  the  previous 
sections,  it  is  concluded  from  these  results,  that  the  early  driver  gas  contamination  is  successfully 
prevented  by  the  sleeve. 


5.1.2  Low  Enthalpy  Operating  Condition  XIII 

As  an  additional  example  of  the  calibration  process  of  HEG  operating  conditions,  the  calibration  of 
condition  XIII  is  discussed  here.  Compared  to  the  conditions  considered  in  section  5.1.1.,  this  condition  is 
termed  a  low  enthalpy  condition  (h0  =  3.3  MJ/kg).  For  the  calibration  measurements  the  rake  shown  in 
Figure  10  was  utilised  allowing  to  determine  static  pressure  profiles  in  addition  to  Pitot  pressure  and 
stagnation  point  heat  transfer  profiles.  For  the  nozzle  flow  computations  using  the  TAU  code,  different 
RANS  turbulence  models  were  applied  along  with  thermal  equilibrium  and  non-equilibrium  computations 
to  determine  the  influence  of  different  modelling  assumptions  on  the  computed  free  stream  conditions.  It 
should  be  noted  that  for  condition  XIII,  the  chemical  relaxation  process  is  in  equilibrium  and  no  free 
stream  dissociation  exists.  The  computational  grid  consisting  of  about  20,000  grid  points  and  the  Mach 
number  contours  resulting  from  the  thermal  non-equilibrium  computation  are  shown  in  Figure  36. 


X  [m] 


Figure  36:  CFD  grid  and  Mach  number  contours  for  the  HEG  nozzle  flow  (operating  condition 

XIII) 


Figure  37:  Comparison  of  measured  and  computed  normalised  Pitot  pressure  (left),  static 
pressure  (middle)  and  stagnation  point  heat  flux  distributions  at  the  nozzle  exit  plane;  NEQ: 
thermal  non-equilibrium,  EQ:  thermal  equilibrium,  SA:  Spalart-Allmaras  turbulence  model,  KW: 
Wilcox  k-w  turbulence  model,  NRP:  Nozzle  reservoir  pressure  (HEG  operating  condition  XIII) 
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A  reservoir-pressure  inflow  boundary  condition  is  used  at  the  subsonic  inflow  of  the  nozzle.  The 
thermodynamic  conditions  at  the  inflow  are  computed  using  isentropic  expansion  from  prescribed  nozzle 
stagnation  conditions  using  the  inflow  velocity  vector  which  is  part  of  the  CFD  solution.  The  nozzle 
supply  conditions  for  the  considered  HEG  operating  conditions  are  given  in  Table  1.  The  chemical  non¬ 
equilibrium  5  species  and  17  reactions  rate  set  for  air  proposed  by  Gupta  ([26])  was  applied.  The 
considered  species  are  molecular  and  atomic  nitrogen  and  oxygen  (N2,  02,  N,  O)  and  nitric  oxide  (NO). 
The  CFD  results  are  subsequently  compared  with  Pitot  pressure  and  stagnation  point  heat  flux 
measurements  on  spherical  probes  and  static  pressure  measurements  obtained  with  the  HEG  calibration 
rake.  The  comparison  between  numerical  and  experimental  calibration  data  also  included  the  computation 
of  the  flow  past  the  probes  used  in  HEG  [42]. 

In  Figure  37,  the  comparison  of  measured  and  computed  normalised  Pitot  pressure,  static  pressure  and 
stagnation  point  heat  transfer  profiles  at  the  nozzle  exit  plane  are  shown.  The  data  is  normalised  using  the 
nozzle  reservoir  pressure  and  a  heat  flux  reference  value  derived  from  the  stagnation  point  heat  flux 
measurements  on  the  spherical  calibration  probes.  Regarding  the  Pitot  pressure,  the  computed  data 
resulting  from  computations  assuming  a  thermal  equilibrium  or  thermal  non-equilibrium  nozzle  expansion 
lie  within  the  experimental  scatter  bars.  For  the  present  operating  condition,  the  wall  boundary  layer  is 
assumed  to  be  fully  turbulent  and  the  difference  due  to  the  application  of  different  turbulence  models  is 
negligible.  The  excellent  reproduction  of  the  measured  Pitot  pressure  profile  shape  confirms  the 
assumption  regarding  the  state  of  the  boundary  layer.  The  best  agreement  between  computed  and  average 
measured  Pitot  pressure  profiles  is  obtained  with  the  thermal  equilibrium  assumption. 

The  computed  static  pressure  profiles  reveal  pronounced  deviations  resulting  from  the  application  of 
different  thermal  relaxation  models.  Excellent  agreement  of  the  numerical  and  experimental  results  was 
achieved  using  the  thermal  equilibrium  assumption.  This  result  highlights  the  importance  of  the  static 
pressure  measurements  since  they  help  to  select  the  most  suitable  numerical  model  to  determine  the  free 
stream  pressure  and  consequently  the  free  stream  temperature. 

The  right  plot  of  Figure  37  shows  the  comparison  of  the  numerical  and  experimental  normalised 
stagnation  point  heat  flux  profiles.  Similar  to  the  Pitot  pressure,  the  computed  values  of  this  quantity 
reveal  only  slight  variations  when  the  assumption  of  a  thermal  equilibrium  or  thermal  non-equilibrium 
nozzle  expansion  is  used.  The  obtained  differences  lie  within  the  experimental  scatter  bars. 

5.2  High  Enthalpy  Cylinder  Shock  Layer  Investigations 

When  considering  the  re-entry  trajectory  of  a  space  vehicle  returning  from  low  earth  orbit  into  the 
atmosphere,  the  most  critical  point  concerning  the  heating  loads  on  the  vehicle  is  found  in  the  continuum 
flow  region  in  approximately  70  km  altitude.  In  this  region  of  the  re-entry  path,  the  velocity  of  the  vehicle 
is  approximately  6  km/s  and  the  flow  past  the  vehicle  is  accompanied  by  strong  shock  waves,  leading  to 
high  temperatures  ensuing  dissociation  reactions.  The  fundamental  influence  of  the  thermal  and  chemical 
relaxation  processes  caused  by  these  high  temperature  effects  on  the  external  aerodynamics,  i.e.  the 
pressure  distribution,  flap  efficiency,  shock/shock  and  shock/boundary  layer  interactions  and  on  the 
heating  loads  can  be  investigated  by  looking  at  the  flow  past  basic  generic  flow  configurations  which  are 
especially  designed  in  order  to  focus  on  one  of  these  effects.  Additionally,  these  studies  are  well  suited  to 
validate  the  ground  based  facility  performance,  measurement  techniques  and  computational  fluid 
dynamics  (CFD)  codes. 

A  test  campaign  was  performed  in  HEG  to  study  the  relaxation  processes  in  the  shock  layer  of  a  cylinder 
placed  with  its  axis  transverse  to  the  flow  [42].  This  configuration  was  chosen  because  of  the  large  shock 
stand-off  distance  that  permits  optical  measurement  techniques  to  investigate  gas  properties  in  the  shock 
layer.  The  diameter  of  the  cylinder  is  90  mm  and  its  length  amounts  to  380  mm.  The  cylinder  model  is 
equipped  with  17  pressure  transducers  to  measure  surface  pressure  distributions  and  17  thermocouples  to 
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measure  surface  heat  flux  distributions.  These  transducers  were  distributed  along  six  rows  located  in  the 
plane  of  symmetry  at  midspan  location  and  10,  20  and  30  mm  to  the  left  and  the  right  of  the  plane  of 
symmetry.  The  model  was  mounted  on  the  nozzle  centreline  (Figure  38). 


Figure  38:  Cylinder  model  in  the  HEG  test  section  including  the  grid  used  for  the  three- 
dimensional  flow  field  computations  (left)  and  model  dimensions  and  gauge  positions  (right) 


Surface  pressure  and  surface  heat  flux  measurements  and  the  determination  of  two  dimensional  shock 
layer  density  distributions  by  phase  step  holographic  interferometry  were  performed.  The  results  which 
will  be  discussed  here  were  obtained  using  HEG  operating  condition  III  (h0  =  12  MJ/kg,  see  Table  1).  Two 
and  three-dimensional  computations  of  the  flow  past  the  cylinder  model  were  performed  assuming  laminar 
flow,  chemical  nonequilibrium,  thermal  equilibrium  and  a  fully  catalytic  cylinder  wall.  The  computational 
grid  used  for  the  3D  simulations  is  shown  in  the  left  part  of  Figure  38.  It  consists  of  33x81x101  grid 
points.  The  symmetry  plane  of  the  3D  grid  was  also  used  for  the  2D  computations  in  order  to  avoid 
discrepancies  caused  by  different  resolutions. 

Two  dimensional  simulations  were  carried  out  using  the  reaction  rate  models  by  Park  [65],  Dunn  and 
Kang  [18]  and  Gupta  et  al.  [26].  The  different  free  stream  conditions  used  for  this  study  result  from  nozzle 
computations,  each  performed  with  one  of  the  reaction  rate  models  mentioned  above.  The  resulting  shock 
standoff  distances  were  compared  with  the  experimental  values  and  a  summary  of  the  results  is  given  in 
Figure  39. 


Reaction  rate  models  used  for  nozzle  flow 
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Reaction  rate  models  used  for  cylinder  flow 


Figure  39:  Comparison  of  the  computed  shock  stand-off  distance  resulting  from  2D  flow  field 

analyses  with  the  values  obtained  in  HEG 


RTO-EN-AVT-1 86 


8-47 


A  Closely  Coupled  Experimental  &  Numerical  Approach  for  Hypersonic  &  High 
Enthalpy  Flow  Investigations  Utilising  the  HEG  Shock  Tunnel  &  the  DLR  TAU  Code 


rrt'B'f-iX-SjT'x; 


It  can  be  seen  that  using  the  reaction  rate  model  of  Gupta  et  al.  results  in  the  smallest  deviation  between 
computed  and  measured  shock  stand-off  distance.  For  that  reason  this  chemistry  model  was  also  chosen 
for  further  2D  and  3D  cylinder  flow  field  computations.  In  order  to  investigate  the  influence  of  the  flow 
past  the  edges  of  the  cylinder  on  the  flow  in  the  central  part  of  the  model,  3D  flow  field  computations  were 
performed.  A  comparison  between  2D  and  3D  results  for  surface  pressure  and  heat  flux  along  the 
stagnation  line  in  span  wise  direction  is  shown  in  Figure  40.  The  centre  plane  of  the  cylinder  is  located  at 
z/R=0.  Two  different  sets  of  inflow  conditions  were  used  for  the  3D  simulations.  The  first  set  is  referred  to 
as  "parallel"  and  represents  a  parallel  inflow  with  averaged  free  stream  conditions  in  the  plane  of  the 
cylinder  position.  The  second  set  is  referred  to  as  "conical".  Here  the  diverging  flow  resulting  from  the 
conical  nozzle  computation  is  used  to  interpolate  inflow  conditions  on  the  CFD  grid  around  the  cylinder. 
A  comparison  of  the  computed  surface  pressure  and  heat  flux  distributions  in  the  centre  plane  is  given  in 
Figure  41.  It  is  clear  from  these  plots  that  the  surface  data  in  the  centre  plane  resulting  from  3D 
computations  compares  well  with  the  data  obtained  from  the  2D  computation.  The  heat  flux  and  pressure 
distributions  along  the  stagnation  line  of  the  cylinder  obtained  with  the  3D  computations  show  the 
influence  of  the  cylinder  edge  effects  on  the  surface  pressure  and  heat  flux.  From  the  centre  plane  up  to 
about  z/R=2.5,  no  significant  impact  of  the  edge  effects  exists.  It  can  also  be  seen  from  Figure  40  that 
compared  to  the  “parallel”  inflow  condition,  the  utilisation  of  the  “conical”  condition  results  in  only  a 
small  difference  of  the  obtained  surface  pressure  and  heat  flux  in  the  centre  plane. 
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Figure  40:  Comparison  of  surface  pressure  and  heat  flux  along  the  stagnation  line  in  spanwise 
direction  resulting  from  2D  and  3D  computations 
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Figure  41 :  Comparison  of  surface  pressure  and  heat  flux  in  the  center-plane  resulting  from  2D 

and  3D  computations 
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The  comparison  between  computed  and  measured  surface  pressure  and  surface  heat  flux  measurements  for 
HEG  run  627  is  shown  in  Figure  42.  Good  agreement  between  the  2D  predictions  and  the  measurements 
was  obtained  (4%  deviation  for  stagnation  pressure  and  1%  deviation  for  stagnation  heat  flux).  Based  on 
previous  investigations  it  is  obvious  that  the  surface  pressure  is  not  very  sensitive  with  respect  to  the 
physico-chemical  modelling  in  CFD  tools.  Therefore  the  potential  of  this  quantity  for  code  validation  can 
be  regarded  as  low.  For  the  present  flow  and  surface  property  conditions,  this  insensitivity  was  also 
observed  for  the  surface  heat  flux.  However,  shock  stand-off  distance,  density  distribution  in  the  shock 
layer  and  static  free  stream  pressure  are  sensitive  to  variations  in  the  physico-chemical  modelling  and  are 
therefore  well  suited  for  validation  purposes. 
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Figure  42:  Comparison  of  measured  and  computed  surface  pressure  and  heat  flux 


Figure  43:  Sketch  of  the  line  of  sight  phase  shift  distribution  reconstruction  based  on  the 
computed  3D  flow  fields  (left)  and  comparison  of  computed  and  measured  phase  shift 
distribution  (right);  HEG  operating  condition  III  (run  627) 
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It  was  shown  through  the  comparison  of  the  results  of  the  two  and  three  dimensional  computations,  that 
the  flow  in  the  symmetry  plane  of  the  cylinder  can  be  regarded  as  two  dimensional.  However  3D  effects 
become  important  for  optical  line  of  sight  methods.  Therefore,  numerical  phase  shift  distributions  were 
obtained  from  the  CFD  solutions  by  using  a  procedure  as  described  by  [31].  A  computational  ray  tracing 
algorithm  was  implemented  (see  Figure  43)  and  numerical  phase  shifts  were  obtained  for  the  2D 
computations  and  the  3D  computations  with  both  “parallel”  and  “conical”  free  stream  conditions.  The 
comparison  of  the  computed  and  measured  phase  shift  along  the  stagnation  line  is  given  in  Figure  44. 
These  Figures  show  clearly  that  the  result  obtained  from  the  3D  computations  by  ray  tracing  along  a  set  of 
lines  of  sight  differ  from  the  2D  CFD  results.  It  should  be  emphasized  that  these  differences  are  only  due 
to  the  contributions  of  the  outer  flow  regions  past  the  cylinder  edges.  The  use  of  the  “conical”  free  stream 
condition  for  the  3D  computation  results  in  a  shock  stand-off  distance  which  is  reduced  by  about  3%. 


Figure  44:  Comparison  of  computed  phase  shift  distributions  along  the  stagnation  line.  The 

phase  shift  values  are  normalized  by  2n 


Figure  45:  Comparison  of  computed  and  measured  phase  shift  distributions  along  the 
stagnation  streamline  for  HEG  operating  condition  III  [54] 

Figure  45  shows  the  computed  and  measured  phase  shift  distributions  along  the  stagnation  streamline. 
Using  the  result  of  the  3D  computation  with  “conical”  free  stream  condition,  very  good  agreement  was 
obtained  with  the  considered  experimental  phase  shift  distribution  [54].  The  corresponding  two- 
dimensional  phase  shift  distributions  are  plotted  in  the  right  part  of  Figure  43. 
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5.3  Free  Jet  Testing  and  CFD  Analysis  of  the  HyShot  II  Scramjet  Flight  Experiment 

In  the  framework  of  the  European  Commission  co-funded  LAPCAT  (Long  term  Advanced  Propulsion 
Concepts  and  Technologies)  project  [86],  the  methodology  of  a  combined  ground  based  testing  and 
numerical  modelling  analysis  of  supersonic  combustion  flow  paths  was  established.  The  approach  is  based 
on  free  jet  testing  of  complete  supersonic  combustion  ramjet  (scramjet)  configurations  consisting  of 
intake,  combustor  and  nozzle  in  the  HEG  and  CFD  studies  utilising  the  DLR  TAU  code  [43],  [53]). 

To  demonstrate  the  capability  of  this  methodology,  the  HyShot  II  flight  test  configuration  ([66],  [82])  was 
selected.  This  generic  scramjet  consisting  of  intake,  combustor  and  nozzle  was  designed  such  that  it 
contains  all  necessary  components  of  a  scramjet  but  still  generates  an  overall  flow  field  which  is  as  simple 
as  possible  with  a  well  defined  combustor  inflow.  Therefore,  it  was  considered  to  be  a  configuration  well 
suited  for  basic  investigations  and  tool  validation  purposes. 

The  HyShot  II  wind  tunnel  model  installed  in  the  HEG  test  section  is  shown  in  Figure  46.  Apart  from  a 
wider  intake,  the  geometry  of  this  model  is  a  1:1  scale  representation  of  the  fuelled  flow  path  of  the 
HyShot  II  scramjet  flight  test  configuration.  Compared  to  the  flight  configuration,  the  width  of  the  present 
intake  ramp  was  extended  from  100  mm  to  196  mm.  The  reason  for  this  is  to  avoid  detrimental  effects  on 
the  flow  visualization  of  the  combustor  flow  by  the  presence  of  external  shocks  close  to  the  outer  surface 
of  the  glass  windows. 


Figure  46:  HyShot  wind  tunnel  model  installed  in  the  HEG  test  section 


A  cut  along  the  symmetry  plane  of  the  HyShot  II  model  configuration  assembly  is  shown  in  Figure  47. 
The  cross  sectional  shape  of  the  intake  and  the  internal  flow  path  of  the  wind  tunnel  model  is  defined  by 
five  parts,  the  intake  ramp,  the  upper  and  lower  combustion  chamber  walls  and  the  upper  and  lower  thrust 
surfaces.  The  intake  is  a  simple  18°  wedge.  Between  the  intake  and  the  combustion  chamber,  a  bleed 
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channel  is  used  such  that  the  boundary  layer  and  entropy  layer  of  the  intake  ramp  do  not  enter  the 
combustor.  Further,  the  bleed  allows  the  shock  generated  by  the  leading  edge  of  the  cowl  to  pass  outside 
the  combustion  chamber.  The  combustion  chamber  is  parallel  and  the  cross  sectional  dimensions  are 
9.8  x  75  mm.  The  length  of  the  constant  area  combustor  is  308.5  mm.  Fuel  is  injected  into  the  combustion 
chamber  by  a  series  of  four  holes,  with  radius  of  1  mm,  equally  spaced  along  one  line  in  cross  flow 
direction  in  the  lower  combustor  wall,  57.5  mm  downstream  of  its  leading  edge. 

Pressures  on  the  model  were  measured  using  KULITE  pressure  transducers.  Wall  heat  transfer  was 
determined  from  the  measured  surface  temperature  histories.  MedTherm  coaxial  thermocouples  were 
utilized.  On  the  intake  ramp,  5  pressure  gauges  and  4  thermocouples  were  integrated  along  the  symmetry 
line.  The  upper  combustion  chamber  wall  was  instrumented  with  23  pressure  gauges  and  1 1 
thermocouples.  The  pressure  transducers  were  positioned  along  the  symmetry  line  and  the  thermocouples 
along  a  line  9  mm  off  the  symmetry  line.  In  the  lower  combustion  chamber  wall,  31  pressure  transducers 
and  15  thermocouples  were  integrated.  The  pressure  gauges  were  positioned  along  the  symmetry  line  and 
the  thermocouples  downstream  of  a  hydrogen  injector  along  a  line  10  mm  off  the  symmetry  line.  The 
upper  and  lower  thrust  surfaces  were  instrumented  with  5  and  6  pressure  transducers,  respectively.  The 
origin  of  the  coordinate  system  used  to  plot  the  measured  and  computed  combustion  chamber  wall 
pressure  and  wall  heat  flux  data  in  the  subsequent  paragraphs  is  the  leading  edge  of  the  lower  combustor 
wall. 

The  computational  domains  and  the  applied  boundary  conditions  used  to  study  the  intake  and  the 
combustor  /  nozzle  flow  fields  of  the  HyShot  II  configuration  are  sketched  in  Figure  48  in  a  plane  cutting 
through  an  injector.  During  the  short  test  time  in  HEG,  no  significant  heating  of  the  model  walls  occurs. 
Therefore,  the  walls  were  assumed  to  be  isothermal  at  a  temperature  of  300  K  and  fully  catalytic.  The 
HEG  free  stream  (condition  XIII,  see  Table  1)  used  for  the  present  studies  is  related  to  a  flight  trajectory 
point  of  HyShot  II  in  an  altitude  of  28  km.  From  the  flight  data  [63],  an  angle  of  attack  of  a  =  3.6°  was 
derived  and  used  for  the  present  studies. 


Figure  48:  Computational  domains  and  boundary  conditions  used  for  the  intake  and  combustor  / 

nozzle  investigations 


In  a  first  step,  a  three-dimensional  (2,500,000  grid  points)  and  a  two  dimensional  computation  (440,000 
grid  points)  of  the  HyShot  II  fuel  off  combustor  flow  were  performed.  The  three-dimensional 
computational  domain  covered  half  of  the  combustor  duct  width  applying  a  symmetry  boundary  condition 
along  the  centre  line.  The  inflow  conditions  were  derived  from  three-dimensional  intake  computations. 
The  computed  and  measured  normalised  surface  pressure  distributions  obtained  on  the  lower  combustor 
wall  are  compared  in  Figure  49.  The  profiles  resulting  from  the  computations  and  from  the  measurements 
agree  very  well.  The  experimental  data  shown  here  represent  mean  values  obtained  from  multiple  runs.  In 
addition  to  the  computed  surface  pressure  data  obtained  on  the  centreline,  streamwise  pressure  profiles 
along  lines  of  25%,  50%  and  75%  combustor  half  width  are  plotted  in  Figure  49  revealing  only  slight 
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variations  in  cross  flow  direction.  Consequently,  the  influence  of  the  combustor  side  walls  on  the  flow 
properties  in  the  vicinity  of  the  symmetry  plane  can  be  regarded  as  small. 

Based  on  this  result,  the  three-dimensional  computational  domain  shown  in  Figure  50  was  selected  for  the 
fuel  on  combustor  computations.  In  order  to  minimize  the  computational  cost  all  existing  flow  symmetries 
were  exploited.  The  combustor  inflow  conditions  result  from  intake  computation.  At  the  outflow  plane  the 
major  part  of  the  flow  is  supersonic  and  consequently  the  conservative  variables  are  extrapolated. 
Symmetry  boundary  conditions  are  used  at  the  span  wise  boundary  planes.  One  is  located  at  z  =  0  mm 
cutting  through  the  centre  of  a  porthole  injector.  The  second  symmetry  plane  is  located  at  z  =  9.375  mm 
representing  half  the  distance  between  two  injectors.  Therefore,  only  one  eighth  of  the  original  span  width 
of  the  HyShot  II  configuration  needs  to  be  modelled.  At  the  bottom  and  top  wall,  viscous  no-slip  boundary 
conditions  are  applied.  Transition  from  laminar  to  turbulent  boundary  layer  flow  at  the  top  and  bottom 
combustor  wall  is  modelled  using  the  criterion  \Ree  /Afe]  =  200  .  Here,  Re0  =  peue6  /  jue  denotes  the 


Figure  49:  Normalised  computed  and  measured  surface  pressure  distributions  on  the  lower 
HyShot  II  combustor  wall;  fuel  off  condition 


Figure  50:  Computational  domain  and  boundary  conditions  used  for  the  three-dimensional 
combustor  fuel  on  computations  (x,y,z  scaled  with  1 :2:2) 

Reynolds  number  based  on  the  boundary  layer  momentum  thickness  and  Me  the  Mach  number  at  the 
boundary  layer  edge.  The  sonic  hydrogen  injection  was  modelled  by  partially  including  the  injector  in  the 
computational  domain.  A  reservoir  inflow  condition  as  used  for  the  HEG  nozzle  computations  was  applied 
at  the  hydrogen  inflow  boundary.  A  total  hydrogen  pressure  of  4.6  bar  and  a  total  temperature  of  300  K 
were  prescribed  resulting  in  a  fuel  equivalence  ratio  of  (j)  =  0.295  .  The  fuel  equivalence  ratio  $  is  defined 

as  the  ratio  of  the  fuel-to-oxidizer  ratio  to  the  stoichiometric  fuel-to-oxidizer  ratio:  $  =  8  •  mH 2  /  rha  . 
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The  computational  grid  consists  of  about  850,000  grid  points.  Structured  prismatic  sub  layers  were  used  at 
the  viscous  walls.  A  dimensionless  wall  spacing  of  y+  =  0(1)  was  used  to  ensure  sufficient  resolution  for 
the  low  Reynolds  number  turbulence  models  and  for  the  computation  of  the  wall  heat  flux.  The  normalised 
wall  pressure  and  heat  flux  distributions  generated  by  hydrogen  combustion  in  the  HyShot  II  combustor 
on  the  lower  wall  are  plotted  in  Figure  5 1 .  The  fuel  equivalence  ratio  was  chosen  such  that  a  steady  state 
combustor  flow  was  obtained.  According  to  the  design  of  the  HyShot  II  flight  experiment,  the  combustor 
is  expected  to  show  a  steady  state  supersonic  combustion  flow  up  to  an  equivalence  ratio  of  (ft  «  0.3  [66]. 
For  clarity,  the  experimental  data  are  subdivided  into  three  groups  with  different  ranges  of  fuel 
equivalence  ratio,  0.266  <  </>  <  0.288,  0.315  <  </>  <  0.329  and  0.341  <  </>  <  0.351.  The  computed 
normalised  pressure  and  heat  flux  distributions  were  obtained  with  the  DLR  TAU  code  using  the  Spalart- 
Allmaras  and  the  Wilcox  k-oo  turbulence  model.  For  all  equivalence  ratios,  the  measured  wall  pressure 
distributions  are  characterized  by  a  continuous  increase  along  the  combustor  indicating  supersonic 
combustion.  Good  agreement  between  the  measured  and  the  computed  pressure  distributions  is  achieved. 
Further,  the  computed  normalised  pressure  distributions  depend  only  weakly  on  the  applied  turbulence 
model. 


Figure  51 :  Normalised  computed  and  measured  surface  pressure  (left)  and  wall  heat  flux  (right) 
distributions  on  the  lower  HyShot  II  combustor  wall;  fuel  on  conditions;  the  numbers  in 
parenthesis  denote  the  HEG  run  numbers 

Regarding  the  comparison  of  measured  and  computed  wall  heat  flux  distributions,  significant  differences 
occur  between  experimental  and  numerical  results.  While  the  computed  heat  flux  increases  continuously 
downstream  of  the  fuel  injection  location,  the  measured  values  remain  approximately  constant  after  an 
initial  stepwise  increase.  A  possible  reason  for  this  discrepancy  between  measured  and  computed  heat  flux 
values  could  be  the  strong  span  wise  heat  flux  gradients  obtained  downstream  of  the  hydrogen  injectors. 
The  heat  flux  gauges  are  positioned  downstream  of  a  fuel  injector  and  even  a  small  flow  misalignment 
could  cause  the  observed  discrepancy. 

In  general,  the  comparison  of  measured  and  computed  pressure  and  heat  flux  distributions  shows  that  the 
pressure  generated  by  supersonic  combustion  can  be  predicted  numerically  with  satisfactory  accuracy. 
However,  apart  from  the  qualitative  difference  between  the  measured  and  computed  heat  flux  distributions 
on  the  lower  combustor  wall,  the  CFD  data  depends  strongly  on  the  applied  turbulence  model.  The 
computed  peak  values  on  the  lower  combustor  wall  differ  by  a  factor  of  almost  two.  This  result  points  out 
that  further  research  is  required  regarding  the  modelling  of  turbulence  in  scramjet  combustors  in  order  to 
reduce  the  uncertainties  of  numerical  wall  heat  flux  predictions.  This  should  include  the  utilisation  of 
RANS,  DES  and  LES  models  to  compute  the  fuel/air  mixing  as  well  as  combustion  including  detailed 
chemistry.  Corresponding  activities  are  planned  in  the  framework  of  the  currently  running  LAPCAT  II 
project  and  the  experiments  performed  as  part  of  LAPCAT  I  will  serve  as  validation  data  base. 
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5.4  Ground  Based  Testing  of  the  DLR  SHEFEX-I  Flight  Experiment 

The  SHarp  Edge  Flight  Experiment  (SHEFEX)  program  of  the  German  Aerospace  Center  (DLR)  is 
primarily  focused  on  the  investigation  of  the  potential  to  utilise  improved  shapes  for  space  vehicles  by 
considering  sharp  edges  and  facetted  surfaces.  The  goal  is  to  set  up  a  future  sky  based  test  facility  to  gain 
knowledge  of  the  physics  of  hypersonic  flow,  complemented  by  numerical  analysis  and  ground  based 
testing.  Further,  the  series  of  SHEFEX  flight  experiments  offers  an  excellent  test  bed  for  new 
technological  concepts  and  in-flight  experimental  sensors  and  it  is  a  source  of  motivation  for  young 
scientist  and  engineers  providing  an  excellent  school  for  future  space-program  managers.  After  the 
successful  first  SHEFEX  flight  in  October  2005,  a  second  is  scheduled  later  this  year  and  additional  flights 
are  planned  for  2015  ff.  A  detailed  discussion  of  the  present  results  and  future  challenges  of  the  DLR 
SHEFEX  program  is  given  by  Longo  [52]. 


Figure  52:  SHEFEX-I  flight  hardware 


Figure  53:  SHEFEX-I  launch  vehicle 


* 


In  the  context  to  identify  and  evaluate  new  technologies  to  advance  hypersonic  space  transportation,  a 
major  goal  of  the  SHEFEX-I  flight  experiment  was  to  exploit  the  enhancement  of  aerodynamic  vehicle 
performance  by  using  sharp  edge  configurations.  Simultaneously,  the  utilisation  of  facetted  thermal 
protection  systems  (TPS)  can  significantly  reduce  the  manufacturing  and  maintenance  cost  compared  to  a 
conventional  TPS  system  consisting  of  thousands  of  individually  shaped  tiles  [51].  Driven  by  the  facetted 
concept,  two  main  criteria  were  used  to  define  the  aerodynamic  shape  of  the  SHEFEX-I  fore  body.  To 
have  as  many  as  possible  facetted  panels  and  to  represent  as  many  as  possible  configuration  details  of 
space  vehicles,  such  as  concave  and  convex  chamfers,  wedged  compression  corners  to  simulate  control 
surfaces  and  sharp  unswept  leading  edges.  As  a  consequence,  the  SHEFEX-I  configuration  became  non- 
symmetric  as  is  shown  in  Figure  52. 

The  launcher  was  a  two-stage  solid  propellant  sounding  rocket  system  conceived  primarily  for  ballistic 
microgravity  experiments.  The  launch  vehicle  consisted  of  a  Brazilian  S30  motor  as  first  stage  and  an 
Improved  Orion  motor  as  second  stage  (Figure  53).  Between  the  facetted  SHEFEX-I  experiment  and  the 
second  stage,  two  cylindrical  modules  were  positioned  which  housed  the  recovery  system,  the  main 
electronics,  the  data  acquisition  devices,  the  power  supply  and  the  cold  gas  attitude  control  system.  Since 


RTO-EN-AVT-1 86 


8-55 


A  Closely  Coupled  Experimental  &  Numerical  Approach  for  Hypersonic  &  High 
Enthalpy  Flow  Investigations  Utilising  the  HEG  Shock  Tunnel  &  the  DLR  TAU  Code 


e-*a‘r±A-*irtv‘ 


the  facetted  body  has  no  aerodynamic  control  devices,  the  second  stage  remained  attached  to  the  fore  body 
until  the  end  of  the  experiment  to  provide  flight  stability  through  its  fins.  The  asymmetry  of  the  payload 
imposed  severe  trim  conditions  for  the  re-entry,  requiring  new  launcher  fins,  stage-adapters  and  ballast 
[52]. 

The  SHEFEX-I  flight  was  launched  from  the  Norwegian  Andoya  rocket  range  in  October  2005.  The 
mission  profile  of  the  SHEFEX-I  flight  and  a  picture  taken  from  the  onboard  camera  are  given  in  Figure 
54.  After  separation  from  the  first  stage,  the  second  stage  was  ignited  in  approximately  20  km  flight 
altitude.  The  burn-out  of  the  second  stage  was  followed  by  the  separation  of  the  nose  cone  at  an  altitude  of 
100  km.  The  nose  cone  protected  the  payload  during  the  launch  phase  of  the  experiment.  Due  to  the 
asymmetric  shape  of  SHEFEX-I  this  was  also  needed  in  order  to  stabilise  the  complete  configuration 
during  launch.  After  the  apogee  in  approximately  211  km  altitude  was  reached,  the  vehicle  started  its 
decent  towards  the  earth.  Experimental  data  was  gathered  between  60  km  and  20  km  altitude.  The  total 
flight  time  was  500  seconds,  including  45  seconds  of  test  time  [90].  The  reconstruction  of  the  flight 
included  the  post  processing  of  the  radar  data  and  based  on  this  the  determination  of  the  atmospheric 
conditions  during  the  flight.  The  reconstruction  of  the  angle  of  attack  a  and  angle  of  yaw  |3  of  the  spinning 
vehicle  was  performed  with  the  help  of  onboard  sensors.  This  assessment  of  the  free  stream  vector  was 
performed  by  using  the  flight  mechanic  data  obtained  with  the  Digital  Miniature  Attitude  Reference 
System  (DMARS)  platform.  The  vehicle  started  re-entering  the  atmosphere  with  an  oscillation  in  pitch  and 
yaw  which  damped  out  over  time  [55].  The  vehicle  finally  achieved  a  stable  flight  attitude  with  a 
decreasing  precession  around  the  flight  vector.  Details  about  the  flight  analysis  are  given  in  [4]  and  [5].  In 
order  to  gain  additional  confidence  in  the  flight  reconstruction  and  analysis,  the  cross-check  with  data 
obtained  during  the  ground  based  testing  post  flight  analysis  and  accompanying  CFD  investigations  is  of 
large  interest.  With  the  SHEFEX-I  flight  and  the  subsequent  numerical  and  experimental  post  flight 
analysis,  DLR  could  for  the  first  time  close  the  loop  between  the  three  major  disciplines  of 
aerothermodynamic  research  -  CFD,  ground  based  testing  and  flight  -  using  own  data. 


Figure  54:  Image  taken  from  the  onboard  camera  during  the  SHEFEX-I  flight  (left)  and  mission 

profile  (right) 

The  flight  reconstruction  in  HEG  was  performed  for  an  altitude  of  20.3  km.  This  trajectory  point  was 
reached  close  to  the  end  of  the  experimental  phase  where  the  oscillatory  motion  of  the  vehicle  became 
small.  HEG  operating  condition  XXI  (see  Table  1)  was  used  to  match  the  flight  Mach  and  Reynolds 
number  as  closely  as  possible.  From  Figure  55  it  can  be  seen  that  the  Reynolds  number  based  on  the 
vehicle  length  agrees  to  2.5%  and  the  Mach  number  to  6%  with  the  selected  flight  point  data. 

The  HEG  stainless  steel  wind  tunnel  model  of  the  SHEFEX-I  configuration  has  a  total  length  of  280  mm 
(see  right  part  of  Figure  56).  The  scaling  factor  to  the  flight  configuration  is  1:2.964.  It  is  instrumented 
with  a  total  number  of  152  thermocouples  and  pressure  transducers.  The  model  was  installed  in  the  HEG 
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test  section  on  a  support  system  which  allows  to  vary  the  angle  of  attack  (+  12°)  and  angle  of  yaw  (±5°) 
(Figure  56). 

For  each  run  in  HEG,  a  numerical  analysis  was  performed  applying  the  DLR  TAU  code.  In  order  to  judge 
the  state  of  the  boundary  layer,  computations  were  performed  assuming  a  laminar  and  a  fully  turbulent 
wall  boundary  layer.  The  Spallart  Allmares  (SAE),  Wilcox  k-co  (WIL),  the  linear  explicit  algebraic  stress 
model  by  Rung  (LEA)  [74]  and  the  shear  stress  transport  model  by  Menter  (MEN)  were  applied.  The 
abbreviations  in  parenthesis  denote  the  labels  used  in  the  figure  legends  below. 


H  [km] 

21.5  21.0  20.5  20.0  19.5  19.0 


Figure  55:  Flight  Mach  number  and  Reynolds  number  compared  to  the  corresponding  values  of 
HEG  operating  condition  XXI;  the  angle  of  attack,  yaw  angle  and  altitude  of  the  selected  flight 

trajectory  point  are  also  given 


Figure  56:  Sketch  of  the  SHEFEX-I  model  integration  in  the  HEG  test  section  (left)  and  model 

design  drawings  (right) 

The  comparison  of  measured  and  computed  normalised  surface  pressure  (left)  and  wall  heat  flux  (right) 
distributions  along  the  lower  SHEFEX-I  symmetry  line  (SYMMB)  are  plotted  Figure  57.  The  pressure 
data  is  normalised  by  the  HEG  reservoir  pressure  p0  and  the  surface  heat  flux  data  is  normalised  with  a 
reference  value  measured  in  the  stagnation  point  of  a  spherical  probe  installed  in  the  HEG  test  section.  The 
diameter  of  this  sphere  is  20  mm. 
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Figure  57:  Comparison  of  measured  and  computed  normalised  surface  pressure  (left)  and  wall 
heat  flux  (right)  distributions  along  the  lower  SHEFEX-I  symmetry  line 


Figure  58:  Comparison  of  measured  and  computed  normalised  surface  pressure  distributions 
along  the  lower  SHEFEX-I  symmetry  line  for  varying  angle  of  attack 

The  error  bars  given  in  the  plots  incorporate  the  measurement  error  as  well  as  the  fluctuation  of  the  signal 
during  the  test  time.  The  measured  and  computed  normalised  surface  pressure  distributions  assuming  a 
fully  turbulent  flow  agree  well.  The  major  difference  between  the  surface  pressure  distributions  resulting 
from  the  laminar  and  turbulent  flow  field  computations  is  a  separated  flow  region  at  the  ramp  hinge  line 
(x=182mm)  occurring  in  the  laminar  flow  solution.  This  can  not  be  identified  in  the  measured  data 
indicating  that  the  boundary  layer  on  the  SHEFEX-I  wind  tunnel  model  is  turbulent.  While  the  difference 
between  the  pressure  distributions  obtained  with  different  turbulence  models  is  small,  the  predicted  heat 
flux  distributions  resulting  from  the  Spallart  Allmares  and  the  Wilcox  k-co,  Rung  and  Menter  SST  model 
show  distinct  differences  (Figure  57).  In  general,  the  best  agreement  with  the  measured  data  was  obtained 
with  the  Spallart  Allmares  turbulence  model.  The  heat  flux  predicted  by  assuming  laminar  boundary  layer 
flow  is  significantly  lower  and  the  distribution  also  differs  qualitatively  due  to  the  flow  separation  at  the 
ramp  hinge  line.  In  Figure  58,  the  measured  and  compute  normalised  wall  pressure  distributions  along  the 
lower  symmetry  line  are  compared  for  varying  angle  of  attack.  Again  good  agreement  was  obtained 
between  the  measured  and  computed  data  and  even  at  the  highest  angle  of  attack  no  flow  separation  at  the 
ramp  hinge  line  was  detected. 
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Figure  59:  Comparison  of  the  HEG  post  flight  analysis  normalised  surface  pressure  data  with 

the  SHEFEX-I  flight  data  base  (a=0°) 

From  Figure  59  it  can  be  seen  that  very  good  agreement  between  the  normalised  surface  pressure 
distribution  along  the  lower  symmetry  line  resulting  from  the  HEG  post  flight  analysis  and  the  SHEFEX-I 
flight  data  base  was  obtained. 
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